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TECHNOMETRICS Fesruary, 1959 


A LETTER FROM THE PRESIDENTS 
OF THE ASQC AND THE ASA 


Early in 1958, the Presidents of the American Society for Quality Control 
and the American Statistical Association were approached with a request to 
consider joint sponsorship of a new publication in statistics. The basis of this 
approach was a desire to assist in meeting certain needs in the area of statistics 
of members of ASQC’s Chemical Division and of ASA’s Section on Physical 
and Engineering Sciences. Each of these groups had previously indicated interest 
in a publication devoted to statistics in the area of experimentation, but it 
appeared that neither group could successfully take such a step alone. 

The reasonableness of the request was recognized and arrangements were 
made to work out details for a joint sponsorship that would be acceptable to 
both societies. Part of this work included a rather complete investigation of 
the possible overlap of the area desired for the new publication with the areas 
now being served. Each society independently came to the conclusion that this 
area could not be covered adequately by any one of the existing facilities. It 
also appeared that the relatively small amount of material of the type under 
consideration being published in other journals would continue. A careful 
study of the circulation potential for the new publication suggested that a 
minimum of 1300 subscriptions could be obtained in the first year and on the 
basis of estimated expenses that it could be made self-supporting by the third 
year. On this basis, the Boards of Directors of the two societies agreed to advance, 
equally, funds sufficient to finance the undertaking. In the case of ASQC, 
monies were made available by the Chemical Division for this specific purpose. 

We recognize the uniqueness of publication under joint control, and the 
attendant responsibility of each organization. For this reason it is an under- 
taking with which the Boards of Directors of both societies are actively 
concerned. Thus, they have concurred in the establishment of a Management 
Committee for the new publication, this committee to consist of three repre- 
sentatives from each society and a chairman. Such a Committee was formally 
established in November 1958, having the following membership: 


Chairman: Paut 8. OtmsteaD, Bell Telephone Laboratories 


Selected by ASQC: J. Y. McCuiure, Convair 
MAYNARD 8. RENNER, Dewey and Almy Chemical Co. 
Harry L. Weurity, Dow Chemical Co. 


Selected by ASA: Irvine W. Burr, Purdue University 
ALMARIN Puituips, University of Virginia 
Donatp C. Rutey, U.S. Bureau of the Budget 





The responsibilities assigned to the Management Committee included: 


1. Selecting a name for the new publication, 

2. Selecting the Editor and Associate Editors, 

3. Establishing the guiding editorial policy of the new publication, 
4. Supervising the publication arrangements, and 

5. Handling the finances of the new publication. 


Actions with respect to each of these responsibilities are, of course, subject 
to review by the Boards of Directors of the two societies. 

At its initial meeting, at which we were privileged to attend, the Management 
Committee took the following principal actions: 


1. The name, TECHNOMETRICS, was selected for the new publication. 
2. The following terse statement of the objectives of the new publication 
was adopted: 


TECHNOMETRICS is to be a journal of statistics for the physical, 
chemical and engineering sciences. 
For the immediate future, this is interpreted as applying to statistical 
experimentation, evaluation and discovery in these sciences. 
. J. Stuart Hunter, Princeton University, was chosen to be the first Editor 
of TECHNOMETRICS for a three year term. 
. Donald C. Riley was chosen to be Treasurer for the TECHNOMETRICS 
account, 
. TECHNOMETRICS will be a quarterly, and be published in February, 
May, August and November. 


The Management Committee, with its continuing responsibility to the 
two sponsoring Societies, will no doubt find it appropriate to take a very active 
part in the publication of this new journal, particularly in its early life. We 
are convinced of the merits of TECHNOMETRICS and feel that physicists, 
chemists and engineers will find it of increasing interest in their fields. We 
assure the Editor and the Management Committee the continued support of 
the two societies and express our confidence in the success of this undertaking. 


C. EvGcrEne FisHer, President 
American Society for Quality Control 


Wa tTeER E. Hoaptey, Jr., President (1958) 
American Statistical Association 


Rensis Likert, President (1959) 
American Statistical Association 
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Response Surface Designs 


for Three Factors at Three Levels 


Rosert M. De Baun 


Stamford Research Laboratories 
d4merican Cyanamid Company 


Stamford, Connecticut 


Many experimental designs have been proposed for exploring an unknown response 
function using a second order approximating polynomial model. Second'order experi- 
mental designs are described for the case of three variables in which each variable 
need only be run at three levels. 


1. INTRODUCTION 


In recent years, practical work has verified the utility of the work of G. E. P. 
Box e¢ al in the study of unknown functional relationships between a continuous 
dependent variable and a number of continuous independent variables. The 
problem and the line of attack can be summarized briefly as follows. 

We imagine a dependent variable » with structure as in (1) where y; represents 
the i-th experimental observation of 7. 


¥=nte (1) 


It is usually, although not necessarily, assumed that the e; are normally and 
independently distributed with mean 0 and variance o*. It is further believed 
that 7 is related to m independent variables (factors) in some unknown fashion 
(2). 

n = f(X, , X2,°°+ Xa) (2) 

It is proposed over a limited region of “factor space” to approximate (2) by 
a polynomial] expansion (3). 


k=m 


y, = Bo+ 2, BX; + 2, BuX;X. + coe sb a (3) 
k=1 

Main interest lies in good prediction of the expectation » over the experi- 
mental region and in the “shape” of the resulting response surface (3). Only a 
secondary interest is taken in the values of the individual coefficients of (3). 
As the orientation of the response surface, as well as its shape, are usually un- 
known, some emphasis has been laid upon the relationship of the variance of 

prediction to the place at which prediction might be desired. 
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It is of course true for all problems in experimental design in continuous 
variables that the way in which the space is defined will be different with different 
experimenters. This arises from differences of informed opinion as to what 
ranges are regarded as being reasonable for each variable, and as to what forms 
of the several variables are most relevant. This is an indeterminacy not in any 
way connected with a method of experimental design. Although unfortunate 
choices of range and coding can of course occur, they will tend to be corrected 
in successive cycles of experimentation. (Box 1957). The experimental factors 
are customarily coded to “design’’ factors with zero mean and equal “spread.” 

Having chosen the region desired to explore, it has seemed reasonable to choose 
experimental points which generated “information” in a symmetrical manner. 
In practice this has been interpreted by choosing designs such that the variance 
of prediction was constant on points equally spaced from the center of the design. 
It must be realized, of course, that this criterion is not of the kind that necessarily 
must be satisfied, but rather, that it gives an objective which might be aimed 
for if it can be attained without creating other difficulties. In practice, designs 
in which the generation of information is not seriously asymmetrical will be 
considered satisfactory. Experimental designs for achieving this aim and related 


aims are discussed in papers already published on this subject (Box and Hunter 
1957). 


2. THREE Factor, THREE LeveL DeEsiens 


All possible designs for three factors at three levels each can be regarded as a 
choice of the numbers n, (p = 1(1)27) where n, is the number of replications 
(0, 1, 2, ---) at the p-th point in the 3° factorial design. A particular example 
of such designs is the three factor, non-central composite design. (Box and Wilson 
1951, Box 1954). In the general response surface situation, however, a higher 
degree of symmetry may be desirable. 

The rotatable response surface designs cf Box and Hunter (1957) require 
that at least five levels of each independent variable be taken. In some cases, 
this is difficult or impossible to achieve. Accordingly, a number of three level, 
three factor designs have been examined for their utility in the situation where 
a second order response surface is to be explored. This communication reports 
the results of this study. ; 

The following subsets of points from the 3° factorial array of 27 points may 
be considered. Each subset is itself symmetrical about the origin. 


a) The 3° factorial itself. 

b) The cube (points at (+1, +1, +1)) 

c) The center point (0, 0, 0) 

d) The octahedron (points at (+1, 0, 0; 0, +1, 0; 0, 0, +1)) 

e) The cuboctahedron (points at (+1, +1, 0; +1, 0, +1; 0, +1, +1)) 


It is proposed to examine second order designs built up of several of these 
sets in combination. 
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In particular, the designs investigated are the following: 


a) The 3° factorial itself. 

b) Cube + Octahedron + x center points. 

c) Cube + 2 Octahedra + n center points. 

d) Cuboctahedron + n center points. 

e) Cube + Cuboetahedron + n center points. 

f) Cuboctahedron + Octahedron + n center points. 


CUBE 
OCTAHEDRON 
CUBC STAHEDRON 
CENTERPOINT 


Figure 1 


3. INFORMATION PROPERTIES OF THE DESIGNS 


Each design has a matrix of sums of squares and cross-products of independent 
variables (X’ X) of the form (4) and the inverse matrix (X’ X)~* clearly is of the 
same form. We denote the elements of (X’X) by a, b, --- f, the corresponding 
elements of (X’X)~* by A, B, --- F. Therefore, A, B, D, and E are given by 
the methods shown in Box and Hunter (1957), while C = c™’ and F = f"'. 
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Following the same paper, the experimental designs are compared to each 
other and to several rotatable designs as follows: Define the weighting function 
I(y) of the regression, given by (5). 


nw = way =[8(¥ 2) 2] (5 


where x is a vector describing the moments of the point in factor space at which 
prediction is desired. In the rotatable designs, I(y) is a function only of distance 
p from the design center (0, 0, 0), no matter what direction is taken (6). In 
this case then, the information distribution can be considered as uni-dimensional, 
and designs compared on a basis of the weight at (0, 0, 0) and its changes along 
the design radii 


p = (Xi + X;+ Xp}. (6) 


The three-level designs, however, have been compared as follows: 

1. The information per point at p = 0, i.e., A~*N~* = I(0, 0, 0)/N. 

2. The distribution of I(y)/N on p along the following radii; the face ray, 
i.e., the ray from the center of the design cube to the center of a cube face 
(X% = p’, X} = Xi = 0); the edge ray, i.e., the ray out to the midpoint 
of a cube edge (Xi = X} = p’/2, Xi = 0) and the corner ray, i.e., that 
out to a cube corner (X? = X? = X3 = p’/3). In this way, we can compare 
the relative non-uniformity of information in these response surfaces. 
The weighting function is given by (7). 


Ig) = (A+ (2B +0) > X?4+ D> X!14+ @E4+F) > Xx)" (1) 


From (6), J(y) will be independent of the separate values of X; for the same 
value of p’? when 


2E+ F = 2D (Rotability criterion) (8) 
I(y) = [A + (2B + C)p? + Dp'y" (8a) 
Accordingly, none of the designs listed are “second-order rotatable.” Table I 
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compares the weight per experimental point along each ray at distances p = 
0, 0.5, 1.0, 1.5 and 2.0 for the several designs mentioned above. We define 
“efficiency” as the value of the prediction (weight, or information), considering 
the total number of experimental points. 

Where the object is to fit a second degree equation, this shows that on a per 
point basis, the 3° factorial is by no means the most efficient of those considered, 
nor is it the most efficient on an absolute basis. It is excelled everywhere by the 
cube plus two octahedra designs with 22 and 24 points respectively. Even its 
virtue of orthogonality is shared with the cuboctahedron design (NV = 16), 
which is almost everywhere a more efficient predictor on a per point basis. 

Of the designs studied, all except the cuboctahedron (VN = 13, 14) and the 
cube plus cuboctahedron (N = 22, 24) are more efficient predictors than the 3° 
factorial. Of the designs containing less than 20 points, the cuboctahedron 
(N = 16) and the cube plus octahedron (NV = 16) are comparable. The first 
design provides orthogonality, 3 d.f. for “true” error and an almost rotatable 
information surface. The second gives more information overall, but in a less 
uniform way, and provides only 1 d.f. for an estimate of true error. The quad- 
ratic coefficients are not mutually orthogonal. 

Of the designs containing 20 or more points, the cube plus cuboctehedron 
(N = 22, 24) provides considerable information in an almost uniform pattern. 
The cube plus two octahedra (N = 22, 24) is a very efficient predictor generally, 
but one whose information distribution is rather non-uniform. It can be seen 
that a response surface situation in three factors, where the experimenter is 
restricted to three levels of each factor, still permits a choice of tolerably efficient 
experimental designs. Such a situation is illustrated in what follows. 


4. ILLUSTRATION—CONFOUNDING IN A THREE-LEVEL DESIGN 


An experimenter wished to examine three continuous factors, but was restricted 
to three levels in each factor. Two other factors (two levels, each qualitative) 
were also under study. Finally there was a sixth factor, whose influence it was 
felt would be only linear. The number of levels of this sixth factor was at the 
choice of the experimenter. Although it was felt that the joint effect of all the 
factors was important, this was a situation where it was also felt desirable to 
consider the effects of the variables separately. This placed a premium on orthog- 
onality. The design selected is shown in Table II. The basic experimental 
design consists of a cuboctahedron plus four center points (De Baun, 1958). 
This design in X, , X, and X; is replicated at each level of the two level factor X,. 
Within runs 1-16, B, is simultaneously partly confounded with B,, , B,; and B,;. 
The reverse confounding pattern is used in runs 17-32. Similarly, within runs 
1-16, Bg is confounded with B,, + Bop + B33 , while the reverse confounding 
pattern is used in the last 16 runs. (De Baun, 1956, Yates, 1937). 

The covariances of the several regression coefficients are all zero except for 
those of B, with B,, , By, and B,;. Analysis of the data of Table II gave the 
following regression coefficients. 
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bo = 43.3625 = —5.7625 bi. = 3.3625 
by 5.3750 0.0250 bis = —1.4375 
be 2.9562 = —2.1250 b., = 0.9000 
b, = 5.8438 


Taste II 
Experimental Design and Data For Response Surface Experiment 
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The orthogonality properties of the design permit us to examine the effect of 
factors X,,X, and X, separately from the response surface. The 95% confidence 
intervals on the coefficients B, , B, and B, are 5.6812 + 0.93, —0.06 + 0.93 
and 2.0375 + 0.47. 
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5. Attias PROPERTIES OF THE DESIGNS 


The altered “moments” of the various designs will produce different alias 
structures. These are illustrated as follows. In all the designs shown the 
expectation of the linear coefficients (b;) can be given as: 


E(b;) = Bs + aBiss + O8s53 + CBee + OBisn (9) 


The alias structure (Box & Wilson) for some of the designs studied is shown in 
Table III. Some marked differences are apparent. The full factorial design 
- has the most favorable bias structure of all. The rotatable arrangement permits 
the cubic coefficient to bias the linear coefficient strongly. Of the designs with 
Jess than 20 points, the cuboctahedron appears to have the most favorable 
alias structure while of the larger designs the cuboctahedron plus octahedron 
has the most favorable alias structure. 


Tas1e III 
Comparison of Experimental Designs—Alias Structure 
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The Analysis of Life Test Data 


R. L. PLAcKEtTT 


University of Liverpool 


The experimental situation described is one in which an external stimulus is 
applied to a sample of n subjects. After a time ¢ a certain number k of the subjects 
have reacted to the stimulus. The problem is then to determine what proportion 
of the population of subjects remain to react after a time T given that 7 >t. Methods 
for solution are described which use order statistics and examples are given. 


1. INTRODUCTION 


Consider an experiment in which a sample of electric lamps is taken, and 
burned with a standard voltage. The object is to determine what is the proportion 
of lamps which will burn for longer than six months. One possibility is to wait 
six months and then see what is the sample proportion still burning, but an 
answer to the question may be required more quickly. Another possibility may 
be to accelerate the testing process, but there is then the difficulty of relating 
performance in the accelerated -test to performance in the standard test. The 
technique adopted here consists in postulating for lamp life a probability density 
which involves unknown parameters, and estimating them from the lives of 
the first few lamps to burn out, say those which last for two months or less. 
The proportion which lasts more than six months can then be calculated from 
the formula for the cumulative probability together with the estimates which 
have been found. We assume that production is in a state of statistical control, 
so that the parameters of the lamp distribution remain constant. 

Although the problem has been described in terms of electric lamps, it occurs 
far more generally. Other fields where predictions are valuable of the proportions 
surviving after a specified time include: radio and television components; textile 
wear; rubber abrasion; operations on cancer patients; labour turnover; business 
failures; and animal experiments. Data which illustrate the kind of variation 
encountered are presented by Berkson and Gage (1952), Davis (1952), Lane and 
Andrew (1955) and Lomax (1954). Such experiments can be described in general 
terms by saying that an external stimulus is applied to a group of test subjects. 
Observation on their times to react ceases either when a specified number k 
have reacted, or after a specified time 7’, which may vary from one subject to 
another. These two kinds of censoring are theoretically distinct, but the differ- 
ences between the corresponding estimates tend steadily to zero as the sample 
size increases. In the methods to be described, we shall assume that k is the 
specified quantity, even though the experiment is terminated after a prescribed 
time. The estimates we give are always unbiased; their efficiency seems likely 
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to be high in small samples and will certainly be very high in large. Essentially 
the same methods will apply when, as is sometimes necessary, the sample is 
censored at both ends. Some transformation of life, say into its logarithm, may 
be necessary to make the assumptions reasonable. 


2. ORDER STATISTICS 


The techniques in the following sections involve ordered samples and the 
expected values of order statistics so that a few preliminary remarks on these 
topics may be helpful. 

Suppose that we take a random sample of five readings from a normal distri- 
bution whose mean is zero and standard deviation one. The result might be: 


—0.84 0.68 0.11 1.05 0.09 


The corresponding ordered sample is obtained by arranging these numbers in 
order of magnitude, which gives 


—0.84 0.09 0.11 0.68 1.05 


If this process is repeated many times and each member of the ordered sample 
is averaged over the values which it takes, we get 


—1.163 —0.495 0.00 0.495 1.163 


which are known as the expected values of the order statistics for a sample of 
five from the standard normal distribution. Tables of them for sample sizes 


n < 50 are given by Fisher and Yates (1948, Table 20) and are reproduced by 
Pearson and Hartley (1954, Table 28), with more detail added for n < 20. 

For a normal distribution whose mean is » and standard deviation co, the 
expected values of the order statistics in a sample of size five are 


pw — 1.16380, pn — 0.495¢, LM, u + 0.495¢, uw + 1.1630 


Thus, given the whole sample, we can estimate o by dividing the sample range 
by 2 X 1.163 = 2.326, which-is just the quantity d, in standard quality control 
chart notation. ° 

Similar remarks apply to other distributions but then of course the expected 
values of the order statistics will be different. Consider the standard logistic 
distribution, which is defined by 


x = log {p/(1 — p)} 


where p is the probability of a value less than zx. This distribution is very similar 
in shape to the normal over the range 0.05 < p < 0.95, and has the advantage 
that a table of logarithms is sufficient for many practical applications of it. 
For a saxaple of size n arranged in increasing order of magnitude, it can be shown 
that the expected value of the r-th member is 


£ = i/@—r+ 1) + Vea —r +2 t+ --- + 1/6 — 1) 


Lives i 
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where we have taken (r — 1) > (n — r). Thus forn = 5andr = 4, the expected 
value is } + 3 = 0.833, as compared with 0.495 for a normal distribution. By 
symmetry, the expected value forn = 5 andr = 2 is —0.833. 

3. GRAPHICAL METHODS 


Suppose that 10 lamps are used in the experiment, and 7 have burned out 
when a report is called for after 2 months = 1464 hours. The 7 lives in hours are: 


1050 1089 1272 1302 1345 1380 1423 


If we assume that lives are normally distributed with mean y» and standard 
deviation o, then the expected values of the ordered lives are respectively 


p — 1.540, nw — 1.00c, un — 0.66¢, pw — 0.380, 
pw —,0.12¢, uw + 0.12¢, u + 0.38¢, 


Lives in hours 
1200 


-| . oO - 0 ‘S$ 
Expected values of order statistics 


Figure 1. 


Using the tables which have already been mentioned. Thus the points 
(-1.54, 1050) (—1.00, 1089) --- (0.38, 1423) are distributed about a line of 
slope ¢ which cuts the axis of lives at the value u. In this way, u and o can be 
estimated by fitting a line to the points by eye, as in Fig. 1. 

We might on the other hand suppose that the lives L had a logistic distribution 


(L — a)/B = log {p/(1 — p)} 
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and in this case the expected values of the ordered lives would be respectively 
a — 2.838, a — 1.728, a — 1.098, a — 0.628, 
a — 0.208, a + 0.208, a + 0.628 


where, for example, 1.09 = 4 + 3 + $3 + 3% + 4. The points (—2.83, 1050), 
(—1.72, 1089), --+ , (0.62, 1423) now give estimates @ and 8, whence the pro- 
portion beyond any specified life L is given by 1/{1 + e'“~**}. Note that 
a = 0 and 6 = 1 for the standard logistic distribution. 

A least squares regression line suggests itself here as providing a better fit 
and more accurate values of yu, ¢, a and 8, but this procedure would not be correct 
because the ordered values are correlated and the degree of correlation must be 
taken into account when fitting the line. Some attention has been given to the 
problem of choosing values of the abscissa so that a least squares fit yields the 
best linear unbiased estimates e.g. Chernoff and Lieberman (1954). This work 
strikes me as beside the point since the purpose of a graphical method is pre- 
sumably to avoid calculations. 


4. Mretuops BAsep ON END OBSERVATIONS 


The fact that the sample range is a fairly efficient estimate of o in small 
samples suggests that we can estimate » and o with reasonable precision simply 
by equating the end observations to their expected values. Thus 


p — 1.54¢ = 1050 
fp + 0.38¢ = 1423 


The solution of this pair of simultaneous equations is 2 = 1349, ¢ = 194. 

By comparing the estimates of this procedure with the best linear unbiased 
estimates, we find the efficiency of f is 89.2% and that of ¢ is 89.7%. For smaller 
samples, the efficiencies are higher: in estimating » and o from the lowest 5 
observations in a sample of 8, the efficiency of i is 95.9% while that of ¢ is 93.3%. 
Correspondingly, these efficiencies steadily decrease as the size of the sample is 
increased while the proportion of the sample which we are using for estimation 
purposes remains fixed. The fall in efficiency is likely to be serious only in large 
samples, where the amount of data available justifies the use of more refined 
techniques, presently to be described. 


5. Exact Metruops In SMatt NorMAL SAMPLES 


We recall that the 7 out of 10 lamps which burned out in two months or less 
gave the following lives in hours 


1050 1089 1272 1302 1345 1380 1423 


The mean of the normal distribution from which these lives are supposed to 
have been taken can be estimated by multiplying the ordered lives by 


0.0244, 0.0636, 0.0818, 0.0961, 0.1089, 0.1207, 0.5045 
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respectively, and adding the products, which gives 2 = 1355. Similarly we can 
estimate the standard deviation by multiplying the ordered lives by 


—0.3253 —0.1757 —0.1058 —0.0502 -—0.0007 0.0470 0.6106 


respectively, and adding the products, which gives ¢ = 200. The appropriate 
coefficients for samples up to size nm = 10, and for any degree of censoring of the 
largest observations, have been calculated by Gupta (1952). His tables have been 
extended to eight decimal places and to nm < 15 by Sarhan and Greenberg 
(1956, 1958). These coefficients give the best linear unbiased estimates of mean 
and standard deviation. A further table of Gupta (Table 5) enables us to calculate 
the standard error of ¢ as ¢*/ (0.1167) = 68. 

We now give a method of estimating » and o from the k smallest observations 
in a sample of n which can be used for normal samples with n < 50 and which 
requires only tables of the normal distribution and the expected values of normal 
order statistics. The notation is explained in Fig. 2, where the curve is a standard 


Ficure 2—Standard normal density function. 


normal density function. Here f, is the density at abscissa ¢, , and p, the area 
to the left of the ordinate; f, is the density at ¢, and q, the area to the right of 
the ordinate. The values ¢, , ¢, , --- , & are the expected values of the first k 
order statistics in a sample of size n, and p, , f; , @ , f, are then all determined. 

The first part of the analysis calls for the calculation of the elements of a 
matrix G according to the formulae in Table 1. A table giving f’/p + tf and 
tf'/p + ?{ — f in terms of ¢ would apply for all sample sizes and eliminate most 
of the labour in this method. 

For example, let us estimate » and o from the 5 smallest observations in a 
sample of 8. The data again refer to lives in hours of electric lamps and they are 


1274 1335 1487 1510 1598 
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TABLE 1 


gu = f7/p. + tf: + 1/2(n + 1) 

gr = 1/(n + 1) 

gue = f2/qe — tfc + 1/2(n + 1) 

ga = hfi/m + tf — fi + t/(n + 1) 
Jer = 2t,/(n + 1) 

gee = tefe/qe — tife + fe + e/(n + 1) 


Here ¢; = —1.42, f; = 0.146, p, = 0.078; and t,; = 0.15, f; = 0.394, g,; = 0.440. 
The remainder of the calculations are set out in Table 2. The estimates f@ and 
é are obtained by solving the equations 


0.8042 — 0.284¢ = 1194 
and 
—0.2692 + 0.848¢ — 237 


The solution is conveniently expressed in the form 


p = 1.402 (1194) + 0.445 (—237) = 1563 
é = 0.469 (1194) + 1.329 (—237) = 217 


Using Gupta’s coefficients, the estimates are 1561 and 216 respectively. 

The calculations also provide standard errors for @ and ¢. In fact, the standard 
error of f is ¢~/(1.402/8) = 91, as compared with 87 using Gupta’s Table 5. 

Supposing therefore that the two experiments which we have analysed are 
carried out on batches of lamps which differ in part of the manufacturing process, 
the corresponding difference in mean life is 1563 — 1355 = 208 hours with a 
standard error of 68° + 917 = 114, suggesting some improvement in per- 
formance. No exact significance test is available here, but the ratio of difference 
to standard error can be taken as approximately normal with zero mean and 
unit variance when no population difference exists. 


TABLE 2 


ts ; i Yi 


—1.42 De nits 

—0.85 ;® Grits 

—0.47 Z NiYi 

—0.15 ) 92iyi 
0.15 


iu it gq 


Totals 
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6. Exact Meruops In SMA.Lt Logistic SAMPLES 


We now give a method of estimating a and 6 from the k smallest observations 
in a sample of n from the logistic distribution 


(L — a)/B = log {p/(1 — p)} 


It can be used for samples of any size but a slightly different method is preferable 
in large samples, as we indicate in Section 7. Only tables of natural logarithms 
are required. The notation in Fig. 3 is the same as in Fig. 2; the only difference 
is that the curve is now a standard logistic density function (i.e. with a = 0 
and 6 = 1). Thus f, is the density at abscissa ¢, , and p, the area to the left of 
the ordinate; f, is the density at ¢, and q, the area to the right of the ordinate. 
The values ¢, , t2 , --- , & are the expected values of the first k order statistics 
in a sample of size n, and are given by the formula 


t= W/a-—-r+)+i/m—r+2)+-:-+1/¢-) 


for (r — 1) > (n — r). Other values are derived by symmetry. 
The coefficients in G are displayed in Table 3, where h denotes 1/(n + 1). 
Again, the preparation of a table would greatly assist the computation of G. 


Fiaure 3—Standard logistic density function. 


TABLE 3 





gu = pill — ~r)(m + A) 

fir = 2p(1 — pr)h 

ge = pel — peXqe + h) 

gu = pi(l — mip — 1) — $h{1 — 2p, — 2hpi(l — pr)} 
gor = —h{l — 2p, — 2t-p(1 — pr)} 

gx = pel — pey(txper + 1) — $h{1 — 2px — A2epr(l — pe)} 
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With the same data as in the previous Section, the corresponding calculations 
are set out in Table 4. The estimates & and 8 satisfy the equations 


0.28353@ — 0.097298 = 430.2 
—0.08097@ + 0.682748 = —40.7 


TABLE 4 


0.01170 
0.03420 
0.04783 
0.05470 


D nits = 
ae guts = 0.68274 


and the solution takes the form 
& = 3.6766 (430.2) + 0.5239 (—40.7) = 1560 


B = 0.4360 (430.2) + 1.5268 (—40.7) = 125 


Standard errors are again provided by this expression for the estimates, and 
the standard errors of & and 6 are respectively 


6/(8.6766/8) and $+/(1.5268/8). 


The parameters a and » both represent the distribution mean, which explains 
why 4 and f are almost identical, but the parameters 6 and oa are not directly 
comparable, except to say that both are measures of scale. It is interesting to 
compare the estimates for proportions of lives greater than a given time which 
are provided by normal and logistic assumptions respectively. The proportion 
of lives beyond 1800 hours obtained by supposing the data are normally distri- 
buted is 0.137, and the proportion on a logistic basis is 0.129. Thus much the 
same prediction is made from the two assumptions. 


7. LARGE NORMAL SAMPLES 


In large samples (n > 50), no tables are available either for the coefficients 
in the linear combinations giving @ and ¢ or for the expected values of order 
statistics. However, a modified version of the general method described in 
Section 5 gives satisfactory results. The modification consists in replacing the 
expected values of order statistics by approximations to them, and it readily 
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extends to situations where the data are grouped, either because observations 
are taken at regular intervals instead of continuously, or because this shortens 
the computing time. 

Suppose that the initial sample consists of 300 lamps, and that the life distri- 
bution of the first 119 is given by the first two columns of Table 5, where e.g. 


TABLE 5 


Central life Frequency Cumulative Average 
(ys) frequency rank 


975 
1025 
1075 
1125 
1175 
1225 
1275 
1325 
1375 
1425 


NRPSaEnNxna ewww 
SSSSESAS awn 
ounmanounounn 
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the entry 975 means that the first group of lives extended over the interval 
from 950 to 1000 hours. These data appear in Gupta (1952). The fourth column 
gives the average rank for each group e.g. the six lives between 1100 and 1150 
have the ranks 8, 9, 10, 11, 12, 13 in ascending order of magnitude and their 
average rank is therefore 10}. The values p; are average ranks divided by 300, 
and the corresponding normal equivalent deviations appear in the column 
headed ¢; . Each value of ¢; thus obtained is taken to represent the average 
expectation of the normal order statistics for the corresponding group. The 
remainder of the analysis follows the pattern already described in Section 5, 
but of course the frequencies in column 2 are now taken into account when 
computing the various sums. We find @ = 1503 and ¢ = 207. These estimates 
compare very well with the values 2 = 1502 and ¢ = 202, obtained by applying 
the principle of maximum likelihood and solving the awkward equations which 
result by the use of a double-entry table specially computed for the purpose. 

The analysis for a logistic distribution combines the relevant features of this 
Section and the preceding one, and any details are therefore superfluous. 

An analysis of grouped frequency data which has been censored is also possible 
by an iterative method due to Hartley (1958) but it seems doubtful whether 
this would be advantageous in the present example where the grouping interval 
is relatively small, and two thirds of the sample are missing. In any case, a 
comparison is difficult because the rapidity with which Hartley’s method con- 


verges depends to some extent on how good are the initial guesses which are 
necessary. 
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8. CoNcLUDING REMARKS 


(a) The foregoing techniques have resulted from an attempt to rationalize 
the methods in this field, which spread far and wide. Independent developments 
have been made on similar lines by Blom (1956, 1958) and Plackett (1958), 
where references to earlier work will be found. 

(b) In practice, attention should be paid to the efficiencies of the censored 
estimates relative to those from the whole sample. It may be more economic 
not to censor but to wait until all lives are completed. 

(c) The numerical methods in Sections 5, 6, and 7 are not very suitable for 
the exponential distribution, which is an adequate description for many instances 
of industrial failure data. However, the relevant estimation and test procedures 
have been studied in detail, notably by Epstein and his co-workers e.g. Epstein 
and Sobel (1953). 

(d) Similar methods apply to censoring in the middle of a distribution, as is 
the case when items fail over the weekend or during an off shift. 

(e) Estimation from incomplete multivariate sample has been looked at by 
several writers, listed by Nicholson (1957) but the techniques which have so 
far resulted hardly seem applicable to the following problem, given in its bivariate 
form. An article is liable to fail from two causes and the joint life distribution is 
specified up to one or more parameters. Experimentation terminates after k 
articles have failed, either from one cause or the other; or at a specified time 7’. 
We want a simple method of estimating the unknowns. 
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1. FREQUENCY DiIsTRIBUTIONS AND RANDOM VARIABLES 


The very existence of quantitative or biometrical work in the Pharmaceutical 
Sciences is in itself evidence of the extent to which during the last sixty or seventy 
years one of the greatest obstacles to exact thought has been effectually overcome. 
The obstacle to which I refer is the existence of variability in the natural world. 
Only one hundred years ago the inhibiting effect of this obstacle can be seen in 
many writers; for us to-day the obstacle does not exist. The familiar concept of 
a frequency distribution as used by Adolphe Quetelet comes to our minds, and 
we recognize that by means of this device, easily represented on paper, varia- 
bility may be accurately specified, and its consequences calculated. 

A close mathematical analogy, though not a logical identity, exists between 
the specification of a frequency distribution, and what in recent years has come 
to be called a Random Variable. If x stand for such a variable and if, for all 
values of P from 0 to 1 we can make probability statements of the form 


Prix < z,) = P, 


where z, is a known function of P, then z is said to be a random variable, and 
knowing z, , we may say it is a random variable of known distribution, for the 
probability of x lying between any two values z, and x, may be equated con- 
ceptually to the frequency in a variable population of objects having some charac- 
teristic measurement which lies between these limits. The probability statements 
would then be equivalent to those appropriate to the knowledge that the subject 
of these statements had been chosen at random from a real population having 
this specification. 

However, the probability statements do not imply the existence of any 
such population in the real world. All that they assert is that the exact nature 
and degree of our uncertainty is just as if we knew it to have been one chosen 
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at random from such a population. The subject of a probability statement if we 
know what we are talking about, is singular and unique; we have some degree 
of uncertainty about its value, and it so happens that we can specify the exact 
nature and extent of our uncertainty by means of the concept of Mathematical 
Probability as developed by the great mathematicians of the 17th century 
Fermat, Pascal, Leibnitz, Bernoulli and their immediate followers. 

The emergence of the notion of mathematical probability in these great 
minds was undoubtedly a major step in the intellectual development of the 
human race. It was unknown to the Greeks, and to the Islamic mathematicians 
of the Middle Ages. It owes its emergence, I suppose, to the high prestige of 
the recreation of gambling among the nobility of France and England, and to 
the existence of a technology advanced enough to supply apparatus of gambling, 
dice, cards, etc. with a precision sufficient to justify the calculations of the 
mathematicians. The type of uncertainty which the new concept was capable of 
specifying with exactitude was that which confronts a gambler in a game played 
fairly, and with perfect apparatus. The probability statements refer to the 
particular throw or to the particular result of shuffling the cards, on which the 
gambler lays his stake. The state of rational uncertainty in which he is placed 
may be equated to that of the different situation which can be imagined in which 
his throw is chosen at random out of an aggregate of throws, or of shufflings, 
which might equally well have occurred, though such aggregates exist only in 
the imagination. 


2. UNCERTAIN KNOWLEDGE OF THE REAL WORLD 


This distinction makes clear the propriety of the application by the great 
Carl Gauss of the Theory of Errors to make probability statements about 
features of the real world, which the observations can ascertain only with 
some uncertainty, as in Astronomy, about the distance of the Sun. Gauss with 
a certain amount of approximation, which later work has removed, showed 
how to calculate a function x, of any value of P from 0 to 1, such that if x stand 
for the unknown distance, the probability statement 


Prix < z,) = P, 


could be asserted for all values of P. He therefore expressed the unknown distance 
x as what we should now call a Random Variable, about which probability 
statements at all levels could be asserted. , 

It is an indication of the state of confusion about even:the most primary 
concepts of mathematical statistics, that more than one modern writer, starting 
perhaps with Neyman about 1935, should have taken exception to this form of 
statement on the ground that, for example, there is only one distance of the Sun, 
that if it is less than x, , the probability should be unity, and if it is greater 
that x, it must be zero. The error, evidently, is to suppose that additional data, 
such as that which provides the exact value of x, can be introduced without 
altering the probability statements which the data actually provide. However, 
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the conclusions of a logical argument must depend on its premises, and it is no 
criticism of the argument that these are different from what they would have 
been had the premises used in the induction been different from what they are. 


3. THE MEANING OF MATHEMATICAL PROBABILITY 


The error has, however, had some troublesome consequences. To obviate 
these and to clarify the matter further, I will set out what I understand to 
be the meaning of the term Mathematical Probability, as used both by the Old 
Masters and by modern scientific practitioners, for it seems to be only in Mathe- 
matical Departments insulated from practical research in the Natural Sciences, 
that confusion and misapprehensions abound. The requirements of a correct 
statement of Mathematical Probability are, I believe, only three: 


(a) A conceptual Reference Set, which may be conceived as a population of 
possibilities of which an exactly known fraction possess some chosen charac- 
teristic. To the extent that this is a meaningful fraction of the whole, the set 
must be measurable, but it need not be measurable in all respects. The set of 
possible throws with a die may be conceived to have exactly 1/6 Aces, and if 
this proportion is the same for throws made on any day of the week, it is un- 
necessary that the specification of the set shall tell us what proportion come 
on Tuesdays. The distribution among days of the week is irrelevant and may 
remain unknown. 

(b) It must be possible to assert that the subject of the probability statement 
belongs to this Set. To the mathematician this may seem trivial, though obviously 
a logical necessity. Tasks of identification, however, belong to the scientist, 
and may require his full attention. We must rely on a Chemist to tell us, if an 
atomic weight has been determined, whether it is that of Potassium or of Rubid- 
ium. This second requirement puts our probability statement into the real 
world. 

(c) No subset can be recognized having a different probability. Such subsets 
must always exist; it is required that no one of them shall be recognizable. 
This is a postulate of ignorance, and therefore unfamiliar to deductive reasoning, 
though characteristic of inductive logic. Let me illustrate it briefly using the 
probability of throwing an Ace. 

(i) The subset of throws made on a Tuesday is easily recognizable, it has, 
however, the same probability as the whole set, and is therefore irrelevant. 

(ii) The subset of throws of odd numbers (1, 3, 5) is obviously relevant, 
since all Aces are odd numbers; this subset is, however, unrecognizable before 
the die is cast. 

(iii) A parapsychologist claiming the gift of precognition might recognize 
a subset in which the Ace is foreseen. If the gift is veridical the proportion of 
Aces in this subset must exceed one in six. In such a case the subset will super- 
sede the original set as the basis of the probability statement. If, however, 
experience shows that the proportion in the subset is just one in six as it would 
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be without precognition, the subset is recognized to be irrelevant, and the 
original set is reinstated. 


Now I submit that these conditions are Necessary and Sufficient; that there 
is nothing lacking in the three conditions I have set out for a valid probability 
statement in the real world; and that there is nothing superfluous. In other 
words, all genuine probability statements fulfil these conditions and any state- 
ment which fulfils these conditions is a genuine statement of mathematical 
probability. 

On this basis we may clear the air somewhat. The concept of Mathematical 
Probability is shown to be unitary. There are not two or more logically distinct 
kinds of probability statements, for if there were they would be susceptible of 
distinct logical specification. Some writers have, however, misunderstood the 
phrase “fiducial probability”, forgetting that it was introduced only to mean 
probability derived by the particular kind of argument to which the term 
fiducial has been applied, and imagining that it signified a special kind of prob- 
ability, for which therefore a special definition would be required. No such 
definition, however, has been proposed. In the minds of these writers we are 
left with a distinction without a difference. 


4. THe FipuctaL ARGUMENT 


I propose to illustrate the fiducial argument in a simple case, which utilizes 
a small but important improvement added in this century to the method of 
Gauss. In 1908 more than 50 years after Gauss’s death, a research chemist, 
W. 8S. Gosset, better known under his pseudonym of “Student”, decided that 
the observations of his laboratory at Guinness’s brewery required for their 
interpretation exact knowledge about the true mean of a normal distribution 
from which only a sample, and perhaps a small sample of observations was 
available. 

His mathematical approach to this problem was to calculate the distribution 
in random samples of the error of the mean divided by the best possible estimate 
of that error. Symbolically, if 1 stand for the true mean, and # for the mean of 
the sample calculated from 

£ = S(x)/N 


where S stands for summation over the sample observed, and if we define s by 
s’ = S(x — Z)’/N(N — 1) 


so that s is the estimated standard error of the mean, then Student determined 
the exact distribution of a quantity ¢ defined so that 


‘(4 — @)/s = t 


The distribution is independent both of the unknown parameter y, and of 
the unknown variance of the population sampled, so that value ¢, can be tabu- 
lated such that for any random sample from any normal distribution 
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Pr(t < t,) = P, 


for all values of P from 0 to 1. 

If, therefore, a chemist has N observations, being independent determinations 
of equal precision of some unknown of importance to his work, he can calculate 
Zand s from his sample, and substitute their numerical values in the expression 
for t so obtaining the probability statement 


Pri{u < (x + st,)} = P. 


The subject of these probability statements is the unknown uy, a property 
of the real world to be determined by observation and experiment like an atomic 
weight, and the series of probability statements about this unique value is 
strictly of the form aimed at by Gauss. 

The improvement on the method of Gauss achieved in the course of 100 
years consists of two elements: (a) Gauss did not know the exact distribution 
of ¢ but took instead as an approximation appropriate for large samples that 
it was normally distributed; (b) In the confusion of the early nineteenth century 
Gauss had indeed introduced into his demonstration the assumption of proba- 
bility a priort supposedly axiomatic as taught by Laplace. He later regretted this 
and in a letter to Bessel recognized its arbitrariness. He was not, however, in a 
position to remedy this logical defect. 

Now, on what conditions do a system of statements such as I have inferred 
by the fiducial argument represent genuine statements of probability in the real 
world? As regards the Reference Set we may recognize that the triad of values 


(x, gz, 8) 


must exist for every sample from every normal population, that for some of 
these samples, but not for all, the inequality 


uw<£+ 8, 


is satisfied, and that whatever may be the true values of the mean and variance 
of the population sampled, the proportion of random samples which satisfy 
the inequality is exactly P. 

The second requirement, that our data really belong to this reference set, 
depends as I emphasized above on the competence of the Scientist and partic- 
ularly on the adequacy of his experimental design. 

The third requirement, which has been the most often ignored, is, however, 
more interesting. Can any subset be recognized within the general set? If a 
distribution @ priori of u were available, this could be used by Bayes’ method 
to recognize a particular subset characterized by 


z + st, 


and to make probability statements appropriate to this subset, different in 
general from those of the fiducial argument. We must therefore specify that if 
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a Bayesian probability a priori is available we shall use the method of Bayes, 
and that the first condition for the applicability of the fiducial argument is 
that no probability a priori of the form needed for Bayes’ theorem shall be 
available. 

Next, it might be thought that some feature calculable from the observations 
could define a relevant subset. However, it has been shown that for Normal 
samples the two statistics Z and s are exhaustive, meaning that the sampling 
distribution of any such function, conditional upon given values of = and s 
shall be completely independent of the true mean and variance of the population 
sampled. There is therefore no relevant and recognizable subset, and all the 
conditions for a genuine statement of probability are satisfied. 

It is sometimes asserted that the fiducial method generally leads to the same 
results as the method of Confidence Intervals. It is difficult to understand 
how this can be so, since it has been firmly laid down that the method of con- 
fidence intervals does not lead to probability statements about the parameters 
of the real world, whereas the fiducial argument exists for this purpose. Moreover, 
the arguments of Neyman and Pearson are deliberately not restricted to cases 
where exhaustive estimation is possible, and this is a definite condition for an 
accurate fiducial argument. It is perhaps only the frequent occurrence of ex- 
haustive estimation in elementary examples which is responsible for the frequent 
coincidence of Confidence with Fiducial limits. The meaning however is always 
different. 


5. BayEes’ THEOREM 


Since the time of Laplace the use by Bayes of probabilities a priori has been 
a stumbling block. For Laplace believed, as Bayes did not, that such probabilities 
were axiomatic a prior. Bayes’ resistance to this notion is shown by two phrases 
in Price’s introductory letter, which explains that in Bayes’ opinion the postulate 
“might not perhaps be looked on by all as reasonable’’, and that therefore he 
had chosen to demonstrate his proposition in another way “rather than to bring 
into his mathematical reasoning anything that might admit dispute.” 

We have these demurrers only at second hand through Price; but if there 
were any doubt we could resolve it by reference to Theorem 8 of Bayes’ own 
text, and see there that the probability statements a priori are not assumed 
axiomatically, but are the result of an auxiliary experiment. 

Bayes’ problem was that presented by Bernoullian or binomial data in which 
a successes and b failures have been recorded out of n trials. His solution took 
the form that if prior to and independently of these observations the unknown 
probability p of success were known as a random variable, so that the probability 
of it falling in the range dp were known to be f(p) dp then the probability of 
this event combined with the outcome of the observations could be expressed as 


es 
aa pa’ f(p) dp, 
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where q stands for 1 — p, and this expression divided by its integral from 0 
to 1, as normalizing factor, would give the probability a posteriori, from which 
probability statements to the effect that p is less than any given value, or lies 
between given limits, would be derived by direct integration. No one in that 
century objected that p could have only one true value, that this must either 
lie between the given limits or not, and that the probability statements were 
therefore meaningless. Such an argument requires a sophistry peculiar to the 
twentieth century! 

How did Bayes obtain his probability a priori? His apparatus was an idealised 
billiard table. He says “I suppose the square table ABCD be so made and levelled 
that if either of the balls O or W be thrown upon it, there shall be the same 
probability that it rests upon any one equal part of the plane as another, and 
that it must necessarily rest somewhere upon it. 

I suppose that the ball W be first thrown and through the point where it rests 
a line shall be drawn’’ parallel to the ends of the table “‘and that afterwards the 
ball O shall be thrown n times’, and that its resting beyond the line should 
constitute a success. 

The auxiliary experiment consists of the first throw, that made with the ball W. 
From this it is inferred that if p be the proportion of the table beyond the line 
drawn, then 


Pr(p < 1) = Di 


for all values of p, , and therefore that the probability that p should lie within 
the range dp must simply be equal to dp, so that f(p) in the analysis above may 
be equated to unity. 


6. A RESTATEMENT OF BAYES’ THEOREM 


The ingenuity of Bayes’ demonstration consists in showing that the single 
probability statement derived from the ball W can be combined with the evidence 
of the ball O, which by itself provides no probability statement about p, to 
give a posteriori an exact expression for p as a random variable. 

There is, however, an element of artificiality about this approach, which 
has perhaps been an obstacle to its understanding. For, it should be noted, after 
the line is drawn we might infer not only the probability a priori of the value p, 
but its actual value, by measuring the distance of the line from the ends of the 
table. 

This artificiality may be easily removed by introducing improved apparatus. 
The billiard table was no doubt sufficiently familiar in the eighteenth century, 
but in the twentieth we can randomize more accurately with a radioactive 
source. 

Let us imagine then a radioactive source emitting a particles at random with 
unknown frequency. This replaces the billiard table. The two balls are replaced 
by detectors of a-particles of two types. The first, which I shall call W, is designed 
to measure the exact time interval between two successive emissions. The first 
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of these starts the clock, the second stops it. The auxiliary experiment with W 
consists in taking a single reading of xz seconds. The probability that it should 
exceed x seconds is represented by P, and it is easily seen that 


-6 
e mz, 


where @ is the unknown rate of emission (in the right direction) supplied by 
the source. Moreover as with Bayes’ ball it is obvious that 


PrP <P) =P,, 


for all values of P, from 0 to 1. 

After the auxiliary experiment therefore, we do not know the value of P, 
but we do know its distribution a priori; what has been actually observed is 
the time interval, z. 

The second detector O, which is to be used n times, is set for a known fixed 
time ¢, and will detect whether in this time interval either no a-particle arrives, 
which will be counted a success, or on the other hand one or more particles are 
detected, which will be a failure. The apparatus is not supposed to be a counter; 
it only distinguishes these two possibilities. The unknown probability of success 
is evidently 


p =e" 


Now the ratio of the times x and é is known, and if 


x/ g= A, 
it follows that 


=P 3 
and 


dP = dp*" dp, 


which is now the probability a.priori of p falling in the range dp. The case 
treated by Bayes is that in which ) is unity, as could easily be realized in practice 
by setting the value ¢ equal to the value of z first observed. 

I suggest, however, that the mild generalization supplied by X is of value in 
reminding mathematicians, who often seem to forget it, that Bayes’ probability 
a priori was obtained, not axiomatically, but by an auxiliary experiment, which 
with the improved apparatus I suggest, could lead to more than one result. 


7. SouRcES OF PROBABILITY a Priori 


Finally, it may perhaps be useful to set out a brief summary of the cases 
which arise in connection with probabilities a priori. 


(a) The commonest case is that in which no probabilities a priori exist as, 
for example, with the chemical determination of an atomic weight. In these 
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cases, when the experiment can be so designed that exhaustive estimation is 
possible, we can infer probability statements a posteriori by the fiducial argument. 

(b) Probabilities a priori may be inherent in the data. This situation is typical 
in Genetics, where before a mouse is born we may be able from its genealogy 
to state in advance the probability of each of its possible genotypes, and to 
combine these probabilities, in the manner of Bayes, with the evidence of the 
frequencies of different types of progeny obtained from test matings. 

(c) It may be possible to follow Bayes’ procedure and to carry out an auxiliary 
experiment which will supply probabilities a priort. The types of observation 
which will do this are indeed the same as will supply by themselves probabilities 
a posteriori by the fiducial argument. 

(d) It has often been imagined that probabilities a priort can be set up axio- 
matically; this has led to many inconsistencies and paradoxes, and I am forced 
to conclude that the process is completely bogus. 


This mid twentieth century is not the first period of grave confusion in the 
teaching of the Theory of Probability. For fifty years from the publication of 
Laplace’s Théorie analytique books in France and England were full of rubbish 
concerning the veracity of witnesses, and the probability of correctness of the 
findings of tribunals. The present confusion seems largely to be a hangover 
from that period from which nineteenth century discussion had largely extricated 
mathematical thought in France and England, but perhaps less completely in 
some more distant countries. 

There is the difference, however, that whereas in the 19th Century error could 
be rife in mathematical departments, without doing greater harm than to confuse 
their students, in our time really important matters such as the standardization 
of drugs, the control of epidemics, and the precision of ballistic missiles are liable 
in the future to be influenced by young men now leaving these departments armed 
with erroneous numerical tables, as well as with confused and obsolete ideas. 
This, in some sort, concerns us all. 
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A Quick, Compact, Two-Sample Test 
To Duckworth’s Specifications* 


Joun W. TUKEY 


Princeton University 


A simple test for comparing two samples is proposed. If one sample contains the 
highest value and the other the lowest value of the two samples combined, the total 
of two counts is taken as the statistic. These counts are of (i) the number of individuals 
in the one sample above all individuals in the other, and (ii) the number of individuals 
in the other below all those in the one. 

5%, 1% and 0.1% critical values are roughly 7, 10 and 13 respectively. Precise 
values for pairs of small and moderate sized samples are given, as are asymptotic 
critical values. Derivations are included. 


I. THe TECHNIQUE IN BRIEF 
1. Introduction 


Since there are a number of procedures for comparing two unpaired samples, 
some which are relatively quick and quite powerful [2, 3, 7], it is natural to 
doubt the usefulness of presenting a new procedure. However, at the Nottingham 
meeting of the Industrial Applications Section of the Royal Statistical Society, 
July 1956, Mr W. E. Duckworth spoke vigorously (in discussion) of the needs 
of certain users for such a procedure which would be very much easier to use 
(and teach) than those so far available. The procedure here presented and dis- 
cussed is an effort to meet his specifications. 

There are procedures which are relatively quick, and highly powerful. To 
merit a place, a new quick procedure must have greater simplicity—and will 
probably have to gain this at the expense of power. Simplicity means practical 
portability—the ability of the statistician to carry the procedure everywhere, 
stored in a very small part of his memory. Such simplicity has to be sought in 
two directions—simplicity in converting the given data, here two samples, 
into a single number (or perhaps two numbers), and simplicity in the table of 
critical values used to assess that number. To be useful, the procedure need be 
not only quick but compact. 

If it is compact, then it can afford somewhat reduced power, for what is 
practically important may be, roughly, the practical power of the procedure 


* Research sponsored by the Office of Naval Research. Tables computed and prepared 
for publication under Contract DA-36-034-ORD 2297, monitored by the Office of Ordnance 
Research. 
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in the sense of Churchill Eisenhart, who has defined practical power as the 
product of the mathematical power by the probability that the procedure will 
be used. A compact procedure may well be used so much more often as to more 
than compensate for its loss of mathematical power. 


2. The most compact form 


Given two groups of measurements, taken under conditions (treatments, 
etc.) A and B, we feel the more confident of our identification of the direction 
of difference the less the groups overlap one another. If one group contains the 
highest value and the other the lowest value, then we may choose (i) to count the 
number of values in the one group exceeding all values in the other, (ii) to count 
the number of values in the other group falling below all those in the one, and 
(iii) to sum these two counts (we require that neither count be zero). If the two 
groups are of roughly the same size, then the critical values of the total count 
are, roughly, 7, 10 and 13, i.e. 7 for a two sided 5% level, 10 for a two sided 1% 
level, and 13 for a two sided 0.1% level. 

(The compactness of the procedure is made clear by the description, complete 
with critical values, in this 110-word paragraph.) 

If the ratio of sizes does not exceed 4 : 3 this procedure will be entirely satis- 
factory. We may extend its usefulness to pairs of samples of less well matched 
size, while preserving semiportability, by subtracting an integer, the correction, 
from the total count before comparing with 7, 10 and 13. Denoting the two 
sample sizes by n and N, where n < N, we define 


0, ifn <N<3+44n/3 
1, if 3 + 4n/3 < N < 2n, 


N-n1n 1 
integer part of Pott , if 2n < N, 


correction = 


Thus for a pair of samples of sizes 7 and 13 the correction is 1, while for a pair 
of samples of sizes 3 and 12, for which (V — n + 1)/n = 3.33 the correction is 3. 
(The fractional part of 3.33 is dropped, leaving only the integer part.) For 
N — n> 9 the critical value 14 is to be used for 0.1% in place of 13. 

The detailed behavior of this form has been examined form < 12,N —n< 11 
(see [9]). The results are summarized in Table 1, both for the correction adopted 
above, and for an alternative correction in which “3 + 4n/3” is replaced by 
hn a = 

A “pocket test’’ of the present sort has rather definite uses. It is for use “as 
a footrule”, ‘‘on the floor’, “in the field’’, etc. Its use is to indicate the weight 
of the evidence roughly. If a delicate and critical decision is to be made, we 
may expect that it will be replaced by some other procedure. Recognition that 
this is the case ought to force us to reconsider our attitude toward the choice 
of critical levels for, procedures based on discrete distributions. If one choice 
of critical value gives a significance level of 3.00%, while the next gives a signifi- 
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TABLE 1 
Behavior of the compact procedure when corrected for inequality of sample sizes. Resulis for 
recommended corrections, n < 12 and N — n < 11. (Results for alternative correction, n < 12 
and N — n <7 given in parentheses.) 


Actual significance levels 


Range Distribution 


Below 3.5% between above 5% 
5% 2.5% to 6.4% 10 80 35 
(2.5% to 5.7%) (26) (61) (2) 
Below 0.7% between above 1.0% 
1% 0.4% to 1.7% 43 44 25 
(0.4% to 1.0%) (59) (15) (0) 


Below 0.07% between above 0.10% 
0.1% 0.06% to 0.20% 12 35 48 


(0.05% to 0.20%) (23) (31) (7) 


Note: Total numbers distributed are not constant because not all nominal 
levels can be reached for certain combinations of n and N. 


cance level of 5.47%, a case which arises with the present procedure, which shall 
we choose for the “‘5% level”? The classical, and conservative, rule is to choose 
the one which comes closest to, but does not exceed, 5%. 

In the case of binomial confidence limits, Hamaker [3] has argued for a ‘‘near- 
est” or “balance of deviations on the average” approach, in the place of this 
“conservative” approach. 

In the present situation it seems rather clearly better to accept 5.47%, instead 
of 3.00% for the nominal 5% probability. The correction set forth above, which 
involves “3 -++- 4n/3’, has been chosen with this in mind, and yields results 
more or less intermediate between those which would be obtained by using 
“conservative” and “nearest” approaches. 


3. Conservative forms. 


For those who want a compact procedure whose significance levels are nearly 
classical in their conservatism, the alternate correction serves quite well for 
N-n< 7. 

When precisely classical significance levels are required, they can be obtained 
from Tables 2A and 2B for the ranges (2A) 0 < N — n < 10 for all n, and 
(2B) 11 < N — n < 20 for N + n < 30. Reliance on these tables obviously 
makes the present procedure much less compact, and hence a poorer competitor 
of existing quick procedures which operate with brief tables of critical levels. 
(In view of the fact that the procedure involves only counting, it may still be 
& quite effective competitor in certain circumstances.) 





JOHN W. TUKEY 


TABLE 2A 


Exact (conservative) two sided critical values for “top count of one sample plus bottom count 
of the other’’, where both counts are > 1, and the sample sizes aren and N withn < N <n + 10. 


N-n n B5%AR/VAGZ|N-—n n S5R/GV/VAM|N—n n 5%/1%/01% 


< 30. 


4 3 8/—/— 8 3 10/13/— 
4-8 7/ 9/13 4-7 8/10/13 4 9/12/— 
9-21 7/10/13 8- 8/10/14 5-6 9/12/15 
7/10/14 7-14 8/11/15 
8/10/14 9/11/— 15-24 8/11/14 
8/11/14 8/10/14 
t}—/— 8/10/14 
7/ 9/— 10/13/— 
7/ 9/13 9/11/— 10/13/16 
7/10/13 8/11/14 9/12/16 
7/10/14 8/10/14 8/12/15 
8/10/14 8/11/15 
3-4 9/12/— 8/11/14 
7/ 9/— 9/12/15 8/10/14 
7/10/— 6-8 8/11/15 
7/10/13 9-17 8/11/14 11/14/— 
7/10/14 18- 8/10/4 10/13/17 
8/10/14. |_-A —_—— 9/13/17 
Note: values for N +n > 30 9/13/16 
7/10/— are extrapolated, using asymp- 9/12/16 
7/10/13 totic results. 8/12/15 
7/10/14 8/11/15 
8/10/4 ‘ 8/11/14 
8/10/14 


bottom count of the other’, 


20 andn + N 


5%/1%/0.1%, for “top count of one sample plus 


4. Analogies 


As a compact, counting-in procedure, the present procedure is part of what 
can now be seen to be a family of actual or needed procedures, namely: 


one sample or two paired samples: the sign test [1] as modified by Duck- 
worth and Wyatt [2] with critical 
difference = 2*/ total number 


two unpaired samples: as given here (with critical values 7, 10 
and 13) 


association between variates: the corner or quadrant sum, [8] with 
critical levels 11, 15(14), and 21(18) 


Exact (conservative) two-sided critical values, 


several samples: nothing yet available (the counting-in 


test [6, 7] is closest, but really requires 
tables) 
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It seems reasonably certain that it would be well to extend this list. The diffi- 
culty lies not in calculating distributions and critical values, but rather in 
selecting simple procedures with reasonably constant critical values. (If a modifica- 
tion of the present procedure could be found which was effective, without correc- 
tion, for samples of markedly different size, it could be extended rather easily 
to a procedure for several samples. If n and N, n < N, are the two sample sizes, 
then the statistic ; 
(end count for sample of n) + 4n/3n + N (end count for sample of N) 

seems to offer promise. It is hoped that someone will investigate this modification 
numerically.) 


5. Example of significance computation 
Given the observations: 
A 147 15.3 161 149 15.1 148 16.7 *17.3 *14.6 15.0 
B 13.9 14.6 14.2 *15.0 14.3 13.8 14.7 14.4 


where we have identified the highest and lowest values in each set with an *, 
we find 


53 A’s above 15.0 (15.3, 16.1, 15.1, 16.7, 17.3 and half of 15.0) 
53 B’s below 14.6 (13.9, 14.2, 14.3, 13.8, 14.4 and half of 14.6) 


for a total count of 5.5 + 5.5 = 11. This reaches the (two-sided) 1% level of 
significance when compared either with the rough 7/10/13 or the precise 
(N — n = 2,n = 8) 7/10/13. 

(As an interim measure at least, the policy of counting tied values one-half 
is recommended.) 


6. Confidence intervals. 


We have discussed the use of the new procedure to test whether “A”’ is sig- 
nificantly less than “B’’, significantly greater, or neither. Often, in fact usually, 
we want to know not only “whether”, but rather “by how much”. We want to 
obtain a confidence interval for the difference between “A” and “B’’. 

The simplest general approach to this problem highlights an assumption, 4 
special case of which we have already made, tacitly, in connection with signifi- 
cance testing. This is the assumption that the distribution of A’s, actual and 
potential, differs only by rigid displacement (“slippage’’) from the distribution 
of B’s, that there is a number 7 so that the distributions of “A — 7” and “B” 
are (substantially) the same. 

If this be so for 7 = m , and we apply the significance test to “‘A — 7’’ and 
“B”, there is but 5% chance of falsely detecting a significant difference. Since 
we do not know 7 , we cannot apply the test for this one value of 7. Instead we 
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must apply it for all values of ». Those values of » for which the test reports 
“not significant” can then be taken as a 95% confidence set for 7, since m will 
be excluded at most 5% of the time. 

For this procedure, as for most others involving a single parameter, the 
confidence set will be an interval. We need then only find out (i) which values 
of 7 are so far to the left that “‘A — 7” is significantly to the right of “B’’, and 
(ii) which values of 7 are so far to the right that “A — 7” is significantly to the 
left of ““B’’. Those 7 not thus disposed of will form the desired confidence interval, 
whose diffidence (= 100% — confidence) will equal the level of significance. 


7. Example of confidence calculation. 


If we wish to deal with the data of Section 5, using 7 as a 5% critical value, 
we might proceed as follows: 


6 = 7 — 1 highest A’s = 17.3, 16.7, 16.1, 15.3, 15.1, 15.0, 
The highest B = 15.0 

(H) <A — B differences = 2.3, 1.7, 1.1, 0.3, 0.1, 0.0, 

6 = 7 — 1 lowest B’s = 13.8, 13.9, 14.2, 14.5, 14.4, 14.6, 
The lowest A = 14.6 

(L) <A — B differences = 0.8, 0.7, 0.4, 0.3, 0.2, 0.0, 


7 highest diffs: 2.3, 1.7, 1.1, 0.8, 0.7, 0.4, 0.3 


If each A were to be reduced by anything less than 0.3, say 0.27, the last 7 
differences would remain positive. Since at least one of these (underscored) 
is from each extreme, “A — 0.27” is significantly higher than “B’’, as will be 
“A — » for any 7 < 0.3. 

We may organize this calculation a little more systematically by arranging 
the two groups (H and L) of differences as follows: 


(either group may be in decreasing order if the other is in increasing order) 


(H) decreasing 2.3 1.7 £3 0.3 0.1 0.0 

(Z) increasing 0.0 0.2 0.3 0.5 0.7 0.8 
lesser of 2: 0.0 02 0.38 O38 0.1 0.0 
highest of these: 0.3 


This format ensures that, if n, be the final value reached, (i) at least 1 difference 
in each group is greater than 7; , (ii) at least 6 + 1 = 7 differences are greater 
than 7, , (iii) is the largest value satisfying (i) and (ii). 

Repeating the construction for B — A differences rather than A — B differ- 
ences, we have (using a conveniently compact format) 
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Lowest 
6 A’s: 14.6, 14.7, 14.8, 14.9, 15.0, 15.1 
1B 138 
Highest 
6 B’s: 15.0, 14.7, 14.6, 14.4, 14.3, 4.2 
1A 17.3 


B-A differences 
(LZ) increasing: —1.3, —1.2, —1.1, —1.0, —0.9, —0.8 


(H) decreasing: —2.3, —2.6, —2.7, —2.9, —3.0, —3.1 
lesser of 2: —2.3, —2.6, —2.7, —2.9, —3.0, —3.1 
highest of these: —2.3. 


thus if 7 is > —(—2.3) = 2.3 we will find ‘A — X” significantly less than “B’” 

Consequently the 95% confidence interval for ‘A — B” obtained by this 
procedure runs from +0.3 to +2.3, that is it contains those 7’s which are not 
<0.3 or >2.3. 


II. Expressions AND CoMPARISONS 
8. General. 


All these derivations are founded on the remark which is fundamental to all 
permutation tests: If the two populations are the same, then (i) the result of 
taking a sample from each and combining their values into a single ordered list 
behaves like an ordered list of a single larger sample from either population, 
and (ii) each order of A’s and B’s in such a list is equally probable. 

We can and will, obtain approximate results for two very large samples 
easily and, with somewhat more effort, exact results for two “‘small’’ samples. 


9. The asymptotic case. 


If the numbers of A’s and B’s in the two samples are in the ratio of p to 4g, 
where p + g = 1, then the chance that any particular value (in the list of the 
two samples ordered together) is an A is p, while the chance that it be a B is q. 
(Note that p = n/(n + N), q = N/(n + N).) 

If both samples are very large, the A-ness or B-ness of one chosen value 
will be very nearly independent of that of any other chosen value. Thus we 
shall have, approximately 


Prob [exactly k highest sre A’s] = p*(1 — p) = p’q 
since the k highest values must be A’s and the next highest value a B. Similarly 
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Prob [at least j lowest are B’s] = qi 


since the j lowest must all be B’s, but the next may be either an A or B. 
How can a total end count of h or more be attained? With A high? Only in 
the following mutually exclusive ways: 


(exactly 1 A high) and (at least h — 1 B’s low) 
(exactly 2 A’s high) and (at least h — 2 B’s low) 


- (exactly h — 2 A’s high) and (at least 2 B’s low) 
(exactly h — 1 A’s high) and (at least 1 B low) 
(at least h A’s high) and (at least 1 B low). 
Hence the probability of a total end count > h with A high is 
(pq)q’* +- (p*g)g’? + +++ + (p**a)a? + (p*"a)a +’) 
= py(qh? + pq? + +++ +p *q +p?) + pg 
Now if p = q = } this becomes 


which must be doubled to obtain the probability of a total end count > h with 
either A or B high. This latter probability is thus h/2* and takes the following 
particular values: 


6 z 9 10 13 14 


h/2’: 3/32 ~ 11% 7/128 ~ 5.5% | 9/512 ~ 1.8% 10/1024 ~ 1.0%} 13/8192 ~ 0.16% 14/16384 ~ 0.09% 


Thus 7, 10, and 14 are at nearly the 5%, 1% and 0.1% levels. 
If, on the other hand, p > gq, the expression in term of p’s and q’s is equal to 


pa'( h-1 + ” (a 4 a) : ro'( 2 oo q* 4 p" = 7m) 
p-4q q q(p — 9) 
h h 
oi‘ion aon 
Pq P—4q 
and the final probability of a total end count of > h with either A or B high is, 
in view of symmetry, 


Pag 

Pg 

which may be shown to reduce to the previous result as p — 3, and is valid for 
p < q by symmetry. 


2pq 
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If we put g = A/(1 + 2), so that X = N/n, 


A N-1 DM d ) 
<-iG47 =i M43 


The latter approximation yields the following asymptotic critical values 


d: (1) 2 3 4 5 9 
5%: (7) 8.1 95 10.6 1.6 14.3 
1%: (10) 12.0 15.0 17.8 20.5 29.6 

0.1% (14) 76 28 Si 61 «28 


The upward drift increases with the extremeness of the significance level, and 
is, perhaps, moderate for 5%. When we allow for use of correction for N ¥ n 
by the rules of Section 2, the approximate 5% asymptotic critical values for 
adjusted counts are 7.1, 7.5, 7.6, 7.6, and 6.3 for \ = 2, 3, 4, 5 and 9. 

Fortunately for the compactness of the test (for pairs of samples of quite 
different size) the situation is more favorable for small and moderate sample 
sizes. 


10. The finite case 


Given n A’s and N B’s, they may be arranged in a variety of orders. The 
total number of ways is specified by “the combinations of N + n things taken 
n at a time” 


(N+")_cNtn_ Waal _ oN tn _(N +n) 


n n n! N! N N 


which is the coefficient of both z* and 2” in the expansion of (1 + x)”*". 
The number of ways for which j > 1, k > 1, 3 + k > h will be shown in 
Section 14 to be given by 


A(N —h,n) — A(N,n — h) 


where A(w, v) is a certain summation of binomial coefficients. Values of A(w, ») 
for small to moderate values of w and v are given in Table 3, which can be easily 
extended by use of the relations (developed in Section 15) 


A(0,v) = 2”"" 
A(w,1) = 1 
A(w + 1,v) = A(w,v) + Aw + 1,0 — 1). 


To determine a two-sided 5%, 1% or 0.1% point, we must calculate, for each 
of a number of values of h, the number of ways in which j >1,k >1,j7 +k 2h 
can be obtained, and compare the results with 


i (N +n) (N+) “ zig (" **) 
40 n |’ 200 “Ft 2000 n /? 


respectively. Table 4 illustrates the process. 
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TABLE 4 
Example of critical value calculation for 
n = 10,N = 20 


A(20 — h, 10) A(20, 10 — h) Difference Comparison 


1,698,160 2,048 1,696,112 
1,097,790 254 1,097,536 


— 751,125 = 


695,860 ; 695,838 
431,910 431,909 
262,144 262,144 
155,382 155,582 


89,946 89,946 
50,643 50,604 
27,824 27,824 


1 /30 
15,022.5 = —— 
ee san 0) 


14,913 _ 14,913 


Resulting critical values: 5% 8 or more 
1% 12 ormore 
0.1% 15 or more 


11. The asymptote as an approximation 


Our tables of exact critical values break off with a suggestion that the asymp- 
totic approximation of Section 9 be used thereafter. Table 5 offers some evidence 
in support of this prescription, showing critical values for N/n = 1, 2, 3, 4 and 
5. The close approximation given by the 5% asymptotic value is clear, as is the 
conservatism of the 1% and 0.1% asymptotic values when used for finite n. 


TABLE 5 
Approach of critical values (5%/1%/0.1%) to their asymptotic value for \ = N/n = constant 


n A=1 A =2 A=3 A=4 A=5 


3 7/10/— 9/11/— 10/13/—  12/15/— 
4 Ippo 8/10/13 9/12/—_—s-:11/15/18 ~=—-:12/17/— 
5 7/ 9/— 8/11/14 9/13/17 10/15/20: 12/18/23 
6 7/9/—- 8/11/14 —-:10/13/17 11/16/21 
7 7/ 9/13 8/11/15 10/14/18 
8 7/ 9/13 8/11/15 10/14/19 
9 7/10/13 8/11/15 10/14/19 
0 7/10/13 8/12/15 10/14/20 
wo  §*/10/14 9*/12/15 10/15/23 11/18/28 12/20/33 


* The next smaller integer barely fails. 
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DETAILS OF DERIVATIONS 
12. Enumeration of ways of obtaining (j, k). 


We noted earlier that there are - -” ways of ordering n A’s among 


N B’s. These ways fall naturally, for our purposes, into 5 ciasses: 
(i) all A’s > all B’s 
(ii) some but not all A’s > all B’s and all A’s >-some but not all B’s 
(iii) highest and lowest either both A’s or both B’s 


(iv) some but not all B’s > all A’s and all B’s > some but not all A’s 
(v) all B’s > all A’s 


Since inverting the order of an arrangement interchanges (i) with (v) and 
(ii) with (iv), and since the arrangements in (iii) cannot lead to total end counts 


which we consider, it will suffice to consider classes (i) and (ii) and then double 
the result. 


Class (i) contains exactly 1 arrangement, with total end count N + n. 
How many arrangements with A end count = j and B end count = k does 
(ii) contain? Any such arrangement must be of the form 


(low) B--- BAC )BA-+- A (high) 


k times j times 


where there are N — j — 1 A’s and n — k — 1 B’s left for the central region. 
These can be arranged in 


pee eee a eee 
n—-k—-1 ro N-j-1 


n—-k-1 


= te GR OL / G—2 ) = ( ire 


~(N-j-D!m—k-1D)! N-j-1 


ways, whereG = N+n—g = N+ 7n — j — k, which must therefore be the 
desired number of arrangements. 

We can thus give the numerical answers in a single table, as is shown in Table 
6, which consists of the customary ‘Pascal triangle’ of binomial coefficients, 
rotated through 45% and with an extra 1 added above the vertex. Thus, for 
example, for VN = 9, n = 7, the number of arrangements with 7 = 3 A’s above 
all B’s and k = 2 B’s below all A’s is found, at coordinates N — 6 = 6 and 
n— k = 5, to be 126. 


13. Enumeration of ways of obtainingg = j +k> h. 


In order to obtain tail probabilities, we need to enumerate the number of ways 
of obtaining j + k > h with j > landk > 1. 


For N = 8 and n = 4, the entries in Table 6 which correspond to possible 
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TABLE 6 


Number of arrangements with j A’s above all B’s and k B’s below all A’s, provided N and 
n are large enough (binomial coefficients specially arranged with an extra 1). 


7 
28 8636 
84 120 
210 330 
462 792 
924 1716 
1716 3432 
3003 6435 12870 
5005 11440 24310 


COON OAR WNH 
eee eee ee 
SOON OAS WN = 


cs 
_ 


arrangements are shown in the upper part of Figure 1. Also shown are the 
entries cut off if we require7 > 1,k >1,g=j+k> 3. 

The remaining entries with the exception of the 1 in the corner fill out a 
pentagonal region, shown in solid lines in the lower part of Figure 1. If we add 


the two triangles A, (composed of values of (“) with v > wu; hence all zero) 


and A, we obtain a parallelogram, over which we will find it easier to sum 
the entries. 

The triangle A’, which we shall use in the argument, represents the increase 
in the pentagonal region if we hold H and N (not h and N) fixed and increase 
n to H or beyond. 

The desired enumeration is now seen to be 


1 + (total for parallelogram) — (total for A) 


not only in this specific case (V = 8, n = 4) but in general. We proceed to find 
these two totals. 

It is natural to sum up the entries in the parallelogram by first summing 
along diagonal strips with g (and G) constant. Such a strip involves values of 
k from 1 to n — 1 and, hence values of n — k — 1 from 0 ton — 2. The corre- 
sponding strip sum amounts (using the result of the previous section, and the 
definition of B(u, v) from Section 14) 
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(75) 4 (774 4..-4+(7 73) = pe@-3,2-9 


The sum over all strips with G = N+n—g<N-—n-—h = His thus 
BO, n — 2) + Bil,n — 2) + --- + B(A — 2,n — 2) 


and by a result of Section 14, this amounts, when the additional 1 is added, to 


BUH — 1,n — 1) 


This, then is the sum over the odd “1” and the parallelogram. 
The sum of the entries over the quadrilateral and the triangle A’ is the result 
of letting m increase to H, namely 


BH — 1,n — 1) = 2" 


Figure 1 


(i) Entries in Table 6 possible for N = 8, n = 4, and losses due to requiring j > 1,k > 1, 
g=jtk>3. 
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whence, using another result of Section 14, 
(sum for A’) = 27" — B(H — 1,n — 1) = BH — 1,H —n-— 1). 

We now recall that H = N+n-—h,sothattH —-N —1l=n-—-h-1, 
whence (sum for pentagonal region) = B(H — 1,” — 1) — B(A — 1), 
n — h — 1) which may be rewritten, with the definitions of Section 14 as 

A(N — h,n) — A(N,n — h) 
since (H — 1) —(n — 1) =H-—-n=N-—h and (A—-1)—(n—-h-1)=N. 
14. Cumulated binomial coefficien's. 


We write (“) for the binomial coefficient, and put 
u u u u 
Bod = ()+(*)+0 2.) + + (9) 


The binomial coefficients are symmetrical, (, ) == (“), so that 


za (e--e (4s jee 


G++ OG t)+- 


= Blu, v) + Blu, u — v — 1) 
which will be convenient later. Since, (") = 0 for u < v, we have 
Buu, v) = Blu,v — 1) = --- = Blu,w), if u<u». 
The fundamental recursion for (“) is 
aaa) 
whence, by repeated use, we find that 
fees Pe eee 


which allows us to complete the marginal sums for the following square array 
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0) er) 
(8) (+) 


eho DOE 
0 0 1 
nd te~ ie Me Reriags + 
where we see that we must have 
Buu + 1,9 +1) = Blu, v) + Bu — 1,0) +--+ + BO,») + 1, 
a result useful in itself, which implies 
Bu + 1,9 + 1) — Bu,v + 1) = Blu, »). 
For earlier convenience, we introduce 
A(w,v) = Bw +» — 1,0 — 1) 
for which 
A(w + 1,v) — A(w,v) = Bw +0, — 1) —- Bwt+v0-—-1,9- 1) 
= Bw +v — 1,0 — 2) 


= Aw+i1,v-—1)), 


AQ,2) = Be —1,0-1)=14+ ("7 ‘) +e 


ad 


A(w, 1) = Bw, 0) = () 


While for w < 0 we have 
- % O0O<wtvc<y, 
A(w, v) = 45 1, wt+v-—0, 
0, w+v<0. 
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Some Statistical Aspects of the 


Economics of Analytical Testing 


O. L. Davies 


Imperial Chemical Industries Limited 


In establishing an economical testing program the costs of routine analytical tests 
must be weighed against losses due to shipping material which is either above or 
below stated specifications. A method for obtaining an optimum sampling and 
testing procedure is described and examples given. 


INTRODUCTION 


The proportion of effort spent on routine analytical testing in a chemical 
factory, and particularly in the manufacture of pharmaceutical products, is 
quite large, and the question of the economics of this testing is a very real one. 
Some reasons for carrying out routine analytical tests are:— 


(i) to check the quality of finished products in relation to the sales specifica- 
tion, 
(ii) to check the quality of raw materials and intermediates for use in manu- 
facturing processes, 
(iii) to estimate the yield of products for control and costing purposes. 


If we did not carry out analytical tests, we would run the risk of serious 
losses; indeed, the purpose of carrying out analytical tests is to avoid these 
losses. The nature of these losses will vary for different situations. For example, 
if a chemical product is released for sale witho. test at an assumed strength 
figure i.e. percent active content, then, because of the variation which must 
exist between batches, some batches will be above strength and (or) some batches 
will be below strength. For those batches released above strength, material 
will be given away unnecessarily and this represents a loss. For those batches 
released below strength, we run the risk of complaints from customers with the 
consequent risk of loss of good-will and loss of business. The latter will result in 
financial loss, although the amount of this loss may be difficult to estimate. 

Even when the strength of the material has been estimated by means of an 
analytical test, these losses still occur, albeit to a lesser extent, because analytical 
tests are subject to error. If the amount of effort spent on testing each batch 
or consignment is increased, the error of the test will be less, and consequently 
the losses in the long run due to releasing material off strength will also be less. 

The situation can be represented graphically as in Figure 1. The horizontal 
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axis represents the amount of effort spent on testing, and the vertical axis the 
losses or costs. For zero effort on testing, the losses due to releasing material off 
strength will be maximal. As the effort on testing increases, so the information 
on each consignment increases and accordingly the resulting losses decrease. 
This is depicted by line Z on Figure 1. 


Relation Between Costs & Analytical Effort 


“ 
oe 
“ 
° 
o 
a 
°o 
“ 
oe 
“ 
a” 
° 
oj 


Effort on Analytical Testing 


Figure 1 


It is also necessary to consider the cost of the effort spent on testing. This 
cost increases with the amount of effort, as depicted by line A: on Figure 1. 
The total costs are represented by curve 7’ which will usually have the shape 
shown on the figure, that is, the curve first of all decreases fairly sharply, passes 
through a minimum, and then increases rather slowly. The most economical 
amount of testing then corresponds to the point of minimum total costs. This 
occurs at the point of balance between the two sorts of costs. If a smaller amount 
of testing effort is used, the losses due to errors in estimation of the quality 
of the material will exceed the saving on testing. If a larger amount of testing 
is carried out the cost of this additional testing will exceed the savings resulting 
from more accurate information on the quality of the material. 

Studies on optimum sampling and testing procedures have been published for 
some years, and a fairly full bibliography is given in the papers by Horsnell 
(Ref 4) and Anscombe (Ref 5). Barnard (Ref 3) has figured prominently in 
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this work. The related problem of the most economical amount of experimenta- 
tion is considered by Yates (Ref 1) and Grundy, Healy and Rees (Ref 2). 

The present paper represents an application and extension of these methods 
to the routine analytical testing of chemical products. - 

In practice there can arise a number of different situations each having its 
own special considerations. I propose to deal with two situations which will 
serve to illustrate the general method of attack. 


First ILLUSTRATION: QUALITY OF FINISHED Propuct 


A product has to conform to a specification with respect to a given quality 
characteristic which for purpose of illustration we will call purity. The specifica- 
tion requires that the purity must exceed a certain prescribed value. Furthermore 
we can measure the purity on a continuous scale with respect to the specification 
limit; and hence give positive values to material above specification and negative 
values -to material below specification. Material below specification can be 
reprocessed to bring it up to specification quality, but a cost is incurred in this 
process. 

For convenience, material above the specification limit will be referred to 
as ‘A’ quality, and material below specification, as ‘B’ quality. 

This classification intoA or B quality is carried out on the result of an analytical 
test. This test is subject to error, and therefore a proportion of the batches of 
material will be wrongly classified. A proportion of A batches will be wrongly 
classified as B quality, and this represents a loss because these batches will be 
reprocessed unnecessarily. A proportion of B quality material will be classified 
as A quality. This suggests a profit, but such a profit could not be maintained 
to any appreciable extent for long because complaints would arise which would 
have to be investigated, and apart from the cost of investigating these complaints 
and compensating those which are justified, there will be a risk of loss of good-will 
and business. Any sort of wrong classification will then result in a loss. 

The loss due to wrongly classifying A material is readily calculated, being 
the cost of reprocessing. The losses which result from selling B material as A 
are not so easy to estimate. The parts which correspond to costs of investigating 
the claims and to compensation could be given a fairly reliable figure, but the 
possible loss of good-will cannot readily be given a cost value. Wherever practic- 
able, attempts should be made to do this, but what we can usually do is to form 
a judgment as to how many such batches can be released without prejudice 
to the sales. For instance it might be possible to conclude that not more than 
one batch in 20 should be wrongly classed as an A, or that we are prepared to 
take a risk of one in 20 of passing as A quality a batch which is, say, 7% below 
the A quality limit. Some such risks are always present, and a person of experi- 
ence who knows the market should be in a position to decide what is a reasonable 
commercial risk. 

To reduce the risk of wrongly classifying a batch, that is to reduce the losses 
involved in such errors, it is necessary to increase the effort on testing. Clearly 
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however it would be uneconomic to increase the testing effort to such an extent 
that the increase in cost of testing exceeds the reduction in the losses due to 
wrong classification. 


OPERATING CHARACTERISTIC 


The risk of wrongly classifying a batch of any given quality can be readily 
calculated when we know the error distribution of the analytical test. Usually, 
the error distribution can be taken as Gaussian with a mean equal to zero and 
a standard deviation o, . The standard error of averages obtained from repeated 
analyses decreases as n increases since 8.E.(g) = o0/ Vn. 

Suppose we carry out a given number n of repeat tests per batch and decide 
to class as A only those batches with a mean test value greater than or equal 
to a given quality X. The probability of classifying a single batch of quality z 
as an A or a B, and accordingly, the probability of wrongly classifying a batch, 
is given by a curve such as that on Figure 2a. This curve is the wellknown 
operating characteristic curve of the test which in this case is given by the 
cumulative Gaussian distribution with mean » = X and a standard error oo/ Vn, 
where op is the standard error of a single analytical test and n the number of 
replicates. 


DISTRIBUTION OF QUALITY 


The distribution of quality from batch to batch can usually be obtained from 
an examination of past records (see Ref 6, p. 73), and the distribution can be 
represented graphically as in Figure 2b. Denote this distribution of the quality 
of incoming batches (measured without error) by f(x) dz. (not infrequently the 
distribution of quality from batch to batch can be taken as Gaussian with a 
mean y and standard deviation o,). Remembering now that positive purity 
values represent material above specification, then the 


Proportion of A batches = [ f(x) dz. 
0 


Now let the probability of classifying as A a batch of quality x be p(z; X, n), 
where n and X are defined above. Then:— 


Proportion of A batches classified as B = 
P= [ (ple; X, mV) ae. () 
0 
If there are any difficulties in evaluating P mathematically, we can use 4 


numerical method similar to that given in Ref 6, p. 77, i.e. we draw ordinates at 
points z; at equal intervals 6, and carry out the summation 


PR » [1 — p(w, ; X, n)If(w.)é 


More accurate quadrature formulae can be used if required. 
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OPERATING CHARACTERISTIC CURVE (a) 
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Optimum ConDITIONS 


The average loss per batch due to wrongly grading A quality batches is simply 


Pc, where c is the cost of reprocessing a batch. The cost of testing per batch 
is bn where b is the cost of a test. 


The losses due to wrongly grading a B are not so readily calculated. However, 
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suppose it was agreed to take a chance a of passing a batch of quality —a; 
this would mean that p(—a; X, n) = a*. The problem then is simply one of 
finding the values of n and X which will minimise the losses Pe + nb subject 
to the restriction 


p(—a;X,n) =a (2) 


A numerical method of doing this is to try out a range of values of n. Substi- 
tuting each n in (2) will enable the corresponding values of X to be determined. 
Given each pair of values X and n, we can evaluate (1). The loss (Pc + nb) is 
then plotted against n and the value of n corresponding to minimum loss is read 
off from the curve. The corresponding X is then obtained from (2). 


NUMERICAL EXAMPLE 


Data (i) The distribution of quality from batch to batch is Gaussian with 

mean y = 0 and standard deviation 7, = 10 units. 

(ii) The distribution of analytical errors is Gaussian’ with mean zero 
and a standard error of a single test o) = 5 units. 

(iii) Cost of reprocessing a batch isc = 500/—. 
Cost of a single test is b = 20/—. 

(iv) A risk of 1 in 20, that is, a = 0.05, is taken of passing a batch with 
quality z = —2.5 units. 


The rule is to carry out n repeat tests, pass as A quality all batches with mean 
result = X and reprocess all other batches. We require to determine the values 
of X and n which will give minimum losses. 

The proportion of batches of A quality which are wrongly classified is given 
by a Normal bivariate integral which can be readily evaluated from tables for 
given values of n and X (see Appendix 1). 

The results of the calculations are as follows:— 





Costs per Batch 
Number Losses due to Acceptance 
of Tests Costsof | wrongly grading Total quality 
=n Testing A quality batches = X 
1 20 103.8 123.8 5.72 
2 40 67.8 107.8 3.31 
3 60 50.1 110.1 2.25 
Ba 80 38.7 118.7 1.61 
6 120 23.9 143.9 0.86 
9 180 11.0 191.0 0.24 


* The principle of fixing one point of the operating characteristic curve and determining 
the other from cost considerations has been applied by Horsnell (Ref 4) to acceptance sampling 
schemes. 
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These results are plotted on Figure 3. The total cost curve has a minimum 
for n = 2 which corresponds to two tests per batch at a cost of 40/—. The 
decision rule is then:— 


Carry out two independent tests on each batch. Average the results, 
and if this gives a value equal to or greater than 3.3, classify the batch 
as an A. Otherwise classify the batch as a B. 


Relation Between Costs & Amount of 
Analytical Testing (First Example). 


200 


Total 


Cost in Shillings 





Number of Tests 


Figure 3 


The cost for n = 3 is only slightly greater than that forn = 2. For n = 3 
the rule requires that batches with a mean test of 2.2 or over should be classified 
as A. Otherwise the batch should be classified as B. 

The operating characteristic curve for n = 2, the distribution of quality from 
batch to batch and the distribution of quality of accepted and rejected batches 
are given in Figure 2. 
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SEQUENTIAL TESTING 


The above method of testing does not exploit fully the possible gains. We 
see from the above example that the standard deviation of batch quality is 
c, = 10 units which is much larger than the standard deviation of a single 
test oo = 5 units. A large proportion of the batches—those which are, say, below 
—10 or above +10—could then be confidently classified with one test, and only 
those batches which are nearer the dividing line between A and B qualities will 
need the larger amount of testing. A sequential testing procedure is therefore 
indicated. 

The simplest procedure of this type is a two-stage test. This would operate 
as follows:— 

Carry out one test on a batch. If the quality is less than X, grade the batch 
as a B. If greater than X; grade the batch as an A, but if the result lies between 
X, and X, carry out two further tests. If the average of all three tests is less 
than X; , grade the batch as a B, but if the average is greater than or equal to 
X; , grade the batch as an A. The quantities X, , X, and X; and if need be the 
numbers of tests at each stage, can be calculated. 

Such a two-stage test will be better than a single stage test. A three-stage 
test will be even better, but the gain in a three or more stage test over a two-stage 
test will probably be quite small. The optimum two, three-stage or sequential 
test can be derived without serious difficulty. Numerical methods for doing 
this will be published separately. 


SECOND EXAMPLE. TESTING OF AN INTERMEDIATE PRODUCT 


This example relates to the amount of testing needed for an intermediate 
product which varies in strength from batch to batch. The intermediate product 
is used in the manufacture of another industrial chemical. The amount of inter- 
mediate to use in a batch of product has to be based on an analytical test, and 
it is assumed that all the other materials are constant from batch to batch. 
If the strength of the intermediate has been overestimated, the batch of final 
product will receive less than its proper charge and the final yield of the batch 
will suffer. If the strength has been underestimated, the batch will receive more 
than its proper charge and there will be some loss due to unbalance with respect 
to the other materials. 

If we plot the cost per lb. of product against the amount of intermediate 
charged, we will get a curve which has a minimum. The point of minimum cost 
represents the optimum process condition with respect to the charge of inter- 
mediate. Any lower or higher charge will result in a loss. It will be assumed 
that the optimum charge is known and has been incorporated in the operating 
instructions of the process. If the optimum is not known, then a simple application 
of Evolutionary Operation should lead one to the optimum, and at the same time 
give the relationship between the cost per lb. and the charge of intermediate 
in the neighbourhood of the optimum, from which the loss per batch resulting 
from any departure from the optimum charge can be readily determined. 
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If there is any error in the estimation of the strength of. the intermediate, 
there will be an error in the charge and this will result in a loss. The degree of 
error can be reduced by increasing the effort (and therefore the cost) of testing, 
and this will result in smaller losses due to errors in the amount of intermediate 
charged. There is an economic limit to the effort spent on testing, because there 
will be a point where any further increase in testing will cost more than the 
amount saved in the cost of the product. What we need is to balance these two 
costs, and this can be done fairly readily. 

The loss curve per batch in the region of the optimum can be represented 
approximately by a quadratic. If x represents the percent error in the charge, 
then the loss will be a constant times 2’, i.e. 


Loss = az” 


Let the standard error of a test be oo , and the cost of a test be c. The accuracy 
increases with the number of repeat tests, the standard error for n repeat tests 
being oo/-Vn and the cost ne. 

When the standard error of estimation of the strength is oo/ Vn, it is shown 


in Appendix 2 that the average loss due to errors in the amounts charged per 
batch will be ac3/n. 


The loss plus the cost of testing will then be (acj/n) + ne. 


This is a very simple expression from which we readily deduce that the mini- 
mum loss is given for n = oo Va/c. 
To take some typical values 


% = 3%, = 10/—, a = 30/— 
then n = 3°/3 & 5. 
It follows therefore that we should spend 50/— on testing per batch. A 
larger effort on testing would cost more than the amount saved through having 


& more accurate charge, and the amount saved by a smaller effort on testing 
would be less than the amount lost in the process. 


The total losses, including cost of testing, for various values of n are:— 


3 4 5 6 7 8 9 10 


120 107.5 104 105 108.6 113.8 120 127 


Although the minimum value is at 5, there is little to choose between the 
values 4, 5, 6 or even 7. But the calculations show that we can lose heavily if 
we rely on one or even two tests. 


UsE or Prior INFORMATION 


The above does not make use of information which may exist on the distri- 
bution of strength from batch to batch. Suppose the prior estimate of the process 
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average is m and the standard deviation from batch to batch is o, . Let the 
testing error standard deviation be oo . When n repeat tests have been carried 
out, the estimate of the strength of the batch based only on these tests is the 
mean, Z, of the replicate tests with standard error oo/‘Vn. But we already have 
an estimate m with standard error ¢, from prior information. These two estimates 
of the process averages are independent and can therefore be combined by 
weighting each by the inverse of its variance. Then 


Estimate of batch strength = (4 m+ a :) / (4, + 4) 
0 


Oy oO} So 


= (oom + nojZ)/(oo + noi) (3) 
The variance of this estimate is o$03/(o05 + no;) (4) 


The method of Appendix 2 can now be applied to obtain the expression for 
the total loss, for all we need do is to substitute the new value of the variance 
for o°:— 


Total loss = ne + aoooi/(o, + no1) 


To find the value of n which gives minimum loss, we differentiate the expres- 
sion, equate to zero and solve for n. This readily gives 


n = (a.Va/e — 09/03) 


The value of n from this expression is, as expected, smaller than the value 
obtained above where the information on the prior distribution was not used. 

m is zero when op = o;V/a/c, i.e. under these circumstances there will be 
no need to test. A less stringent requirement which permits no testing is given 
by equating the total loss expressions for n = 0 and n = 1. However, some 
tests, not necessarily on every batch, will probably have to be carried out to 
detect changes in the average strength and its variation, i.e. in m and o, . This 
is the function of Statistical Quality Control. Naturally, when such tests are 
carried out, the information should be combined with the prior information 
in order to get the best estimate of the strength for each batch. 

The procedure of combining the information on the prior distribution and 
the test results in order to obtain a better estimate, would not, of course, be 
used if there is any reason other than the result of the test, for believing that the 
particular batch is exceptional with respect to the property measured. 

I am indebted to W. Spendley for some useful suggestions, particularly in 
relation to the use of prior information. 


APPENDIX 1 
First Illustration. Quality of Finished Product 


Notation: 
(i) Distribution of quality from batch to batch = f(y) dy. 
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(ii) Distribution of errors due to analytical testing, for fixed quality, are 
assumed to be Gaussian with average zero and standard error ao . 

(iii) Cost per test = b. 

(iv) Cost of reprocessing a batch = c. 

(v) Accepted risk of passing a batch of quality a units less than the specifi- 
cation quality limit = a. 

(vi) The quality is measured with respect to the specification limit taken as 
origin, and all batches with a quality measure = 0 are given a rating of A. 


The decision rule will be of the form:— 


Carry out n repeat tests; pass as A quality all those batches for 
which the mean test result is = X, and reprocess those batches for 
which the mean test result is < X. 
The problem is to find the values of » and X which will give minimum losses. 
For n repeat analytical tests the chance of rejecting a batch of quality equal 
to y is the normal integral. 


pa fr Bp, (—42 de (1) 
2r J-o 


Condition (v) means that (1 — p) = a for y = —a. This gives the numerical 
value of (X + a)-‘Vn/o» from which for any given n, we can find the correspond- 
ing value of X. 


The proportion of batches of quality y is f(y) dy, and therefore the proportion 
of batches wrongly rejected is 


P= [ pi) dy (2) 


where p is given by (1) and the lower limit of the integral is fixed by condition 
(vi) 
*. Total losses = Pc + nb. (3) 


If f(y) is Gaussian with mean m and standard deviation a, , then 


P= Tey [ rE [—(y — m)*/20;] dy 


val” Exp [—(y — m) ‘pail 7 Exp (—2°n/2c5) ax | dy 


~ Qeaoo, 


~ va 5 Exp [—y’/20; — (x — y)’n/2o5] dy dx 


21000; 


This is the Normal bivariate integral taken from (—m/o,) to ~ with respect 


to y and from —- to X/V/(c; + o3/n) with respect to x. The correlation 
coefficient of the surface is o,//o; + o4/n. 


The integral can be readily evaluated from tables for any given value of n 
and its corresponding value of X. 
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Suppose that instead of condition (v), we wish to control the proportion of 
batches passed with quality below zero. 


Proportion of batches passed = 1 — / pf(y) dy = P, 


Proportion of batches passed below specification 


= [0 -pi@ a =P, 


We then have to make P./P; = a. 
Since both P, and P, are Normal bivariate integrals, the value of X for a 
given value of n can be readily determined using the Normal bivariate tables. 


Example 


For the particular numerical example given earlier, we obtain the following 
results:— 


Number of Acceptance Losses due to 
Tests Quality limit Cost of |= wrongly grading Total loss 
=n = X test A batches per batch 


94 114 
55 95 
41 101 
30 110 
19 139 


The optimum number of tests is 2, but there is little to choose between n = 2 
or 3. For the former, the rule would be to pass all batches with a mean quality 
= 2.39, and for the latter, to pass all batches with mean quality = 1.58. 

In this particular example the condition that the proportion of batches of 
B quality passed is equal to 1/20, is less stringent than when the probability 
of passing a batch of quality —2} is made equal to 0.05. The former is probably 
the more useful in practice. 


APPENDIX 2 
Economic amount of testing for intermediate 


The loss per batch of product due to a deviation of x% in the amount charged 
per batch from the optimum is az’, and the standard error per test is c. Then 
the average loss per batch due to testing errors for n replicate tests per batch is 


x / a a 


7 : ° : . . 
— 2; i Oe aoe Sea i a oe et 
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Total loss including cost of testing, is then 


ao’ /n + ne 


where c is the cost of a single test. This total loss is clearly a minimum when 
n= oVa/e. 
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Partial Duplication of Factorial Experiments 


O. Dyxstra, Jr. 
General Foods Research Center 


Tarrytown, N. Y. 


Two level factorial and fractional factorial designs are described in which a subset 
of the treatment combinations are duplicated. The duplicate runs provide an unbiased 
estimate of the experimental error variance and more reliable estimates of the effects. 
Several analysis procedures are described and a numerical example is given. Designs 
for partially duplicated two level factorial and fractional factorial designs are pro- 
posed for as many as eleven factors. 


1. INTRODUCTION 


At the 1957 Convention of the American Society for Quality Control, C. 
Daniel [4] discussed error estimation in factorial experiments, by the device of 
duplicating a portion of the factorial design. It is the purpose of this paper to 
continue this discussion. 


2. DUPLICATION 


As factorial experiments are run in industry, it is often difficult to reset 
temperatures, valve settings, and the like. On the other hand, once having 
attained the experimental conditions prescribed by a design, it is relatively 
easy to obtain a series of samples. The scatter among the resulting measurements 
provides a measure of error, but this error is not the correct one to use in evaluat- 
ing the factor effects. Obtaining duplicates in this manner is essentially the same 
as taking repeated measurements on one sample. However, statements concerning 
the effects of changes in the factors (temperature, pressure, etc.) must be meas- 
ured against the experimenter’s ability to reset and reattain the desired experi- 
mental conditions, i.e., the experimental error, as opposed to the errors of repeated 
observations taken under a particular setting of the factors under study. In 
the type of experiment we envision, the experimental conditions would be 
changed after each run, even though the complete experimental program would 
contain several runs at some of these experimental conditions. In this way we 
obtain what we call ‘‘real duplicates” in the following text. 


3. THe ARGUMENT FOR DUPLICATION IN FACTORIAL EXPERIMENTS 


In factorial experiments where each factorial combination is run once, it is 
customary to determine an “error” estimate by pooling high-order interactions, 
or by means of half-normal plotting, as discussed by Daniel [3]. While pooling 
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the three-factor and higher-order interactions is generally acceptable, we find 
too often that two-factor interactions are also used. As Daniel points out, ‘‘of 
course whenever a two- or three-factor interaction looks big, we quietly snatch 
it from the pool and test its ‘significance’ ’’. In many cases we cannot be certain 
by how much the “error” term is biassed. One solution to this probler is to 
obtain some real duplicates, which provide an unbiassed estimate of the error. 


4. A SIMPLE ILLUSTRATION 
4a. Experimental Plan 


To illustrate partial duplication we use a 2° factorial. The treatment combina- 
tions are (1), a, b, ab, c, ac, be, and abc, using the usual notation. (As a word of 
explanation: The factors are A, B, and C; a factor at its high level is indicated 
by the presence of the corresponding lower case letter. A factor at its low level 
is indicated by the absence of the lower case letter. Thus “‘abc’’ indicates the 
factors A, B, and C at their high levels and “(1)” indicates the factors A, B, 
and C at their low levels.) We will duplicate four of the treatment combinations 
by adding a 2°", that is, a half replicate of the 2° factorial design. In order to 
find the 2°" we confound the ABC interaction, so that we may run the combina- 
tions (1), ab, ac, and be, or a, b, c, and abe. Since there is little to choose between 
these sets, for convenience and simplicity the set is duplicated which includes 
the combination (1), so that J = —ABC is the defining contrast. The complete 
design therefore includes twelve experiments, the factorial combinations a, }, c, 
and abe once and the combinations (1), ab, ac, and be twice. 


4b. Method of Analysis 


From the results of the 2° estimates of the main effects and two-factor inter- 
actions can be obtained, and from the differences between the duplicate runs 
an estimate of the error variance can be obtained. In addition, the results for 
the 2°-* when combined with the results for the 2° provide more reliable estimates 
of the effects. 

The A contrast, represented as (A), is the difference between the sum of all 
results at high A and the sum of all results at low A. While we would usually 
think of the A effect as the average difference between the results at high A 
and the results at low A, i.e., the (A) contrast divided by 4 in the case of a 2’, 
it will be easier to consider the effects as deviations from the overall mean. 
Representing the estimates of the effects as A, B, AB, etc., the averages for the 
levels of A are 7 + A at the high level and 7 — A at the low level. For the 2° 
then, we have A = (A)/8 or (A) = 8A. Similarly for the 2°-* we have (A) = 
4A — 4BC, because of confounding. For the composite design, the 2° and the 
2°", we have (A) = 124 — 4BC. Thus the (A) contrast measures the main 
effect of A and the BC interaction. Similarly, the (BC) contrast_ measures the 
BC interaction and the A main effect. The pair of effects A and BC are given by 
the equations 
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(A) = 124 — 480 
(BC) = —44 + 1280 
The solution is readily found to be 
324 = 3(A) + (BC) 
32BC = (A) + 3(BC) 


Similar equations are found for the pairs of effects B and AC, C and AB, and 7 
and ABC. 

A method of analysis of these data is to compute the sums of squares for 
pairs of degrees of freedom, so that we may use the standard tests of significance. 
The two degree of freedom sum of squares for A and BC is given by 


S.8.for 4 and BC = A(A) + BC(BC) (3) 


(2) 


In a similar manner we may obtain the sums of squares for the other pairs of 
effects. If a mean square is statistically significant, however, we cannot be 
certain which effect is real or whether both effects are real. 

One approach in determining which effects are real, suggested by one of the 
referees, involves a partitioning of the 2 d.f. sum of squares, first assuming 


BC = O and then assuming A = 0. The decomposition of the 2 degrees of freedom 
would be as follows: 


Source df. Sum of Squares 


A and BC 2 A(A) + BC(BC) 
Decomposition 1 


A, assuming BC = 0 (A}*/12 
BC (by subtraction) ((A) + 3(BC)]?/96 


Decomposition 2 


BC, assuming A = 0 (BC)*/12 
A (by subtraction) 1 [3(A) + (BC)I2/96 


The decompositions of the 2 d.f. sums of squares into single degrees of freedom 
should reveal which effects or interactions are real. For the pair 7 and ABC we 
would immediately partition the 2 d.f. into one for the correction factor (c.f. = 
(I)*/12, where (J) is the grand total of the observations) and one for ABC. 
In another approach, given in Section 4d below, the effects are tested sepa- 
rately, but this approach involves ignoring the fact that the effects are correlated. 


4c. Calculation of Contrasts 


Yates [6] gives a means of simultaneously computing all effect-contrasts for 
a 2-series factorial design. The details of the calculations are easily available in 
Bennett and Franklin [1], in Davies [2], or in Kempthorne [5]. 
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It is easily shown that the cuntrasts for a partially duplicated factorial can 
be obtained by using the sums of pairs of observations where they occur and 
otherwise using the single observations. The contrasts which we obtain are 
those which we use in the analysis. 


4d. Tests of Significance 


For a 2” factorial which is half duplicated, the variance of an effect is 3 o°/2"*?, 
if the effect is partially confounded with another effect. Otherwise the variance 
of an effect is o”/3(2""*). Using s’ as the experimental estimate of o” with 2°” 
degrees of freedom, we may test for the statistical significance of an effect F, 
say, by computing 


t = B/VVar (B) (4) 


The effect E is judged to be real if ¢ > ¢, , where ¢, is the value of Student's 
distribution with 2””* d.f. at the a level of significance. 

In performing such tests of significance, we must bear in mind that the effects 
are correlated in pairs, the correlation coefficient being p = }. 


5. AN ALTERNATE METHOD oF ANALYSIS 


It is obvious that we could simply have averaged the duplicate values and 
obtain unbiassed estimates of the effects. For the 2° and 2°-* we had defined 
(A) and (BC) as representing the A and BC contrasts. If, now, we define (A)’ 
and (BC)’ as the A and BC contrasts based on averages of duplicates where 
they appear, it is found that 


(A) = (A)’ + [(A)’ — (BC)’/2, or 2A) = 3(A)’ — (BC)’ 
(BC) = (BC)’ + [(BC)’ — (A)’]/2, or 2(BC) = —(A)’ + 3(BC)’ 
We therefore find 


(5) 


A 


(A)’ = [3(A) + (BC)]/4 = 8A 


“a (6) 
(BCY'= (A) + 3(B0)]/4 = 8BC 


which proves that exactly the same estimates result from either method of 
analysis. Since the estimates are the same, the variances and covariances are 
still the same, and the effects are correlated in pairs. The usual factorial analysis 
based on averages where pairs occur would obviously be in error. 


6. AN ApproacH Usinc Matrix ALGEBRA 


The pairs of equations for the 2° and the 2°~* are derivable directly from a 
regression approach. We define 2° = 8 variables x) through 2; in Table 1. 
The levels of A, B and C are denoted by 2, , x. , and 2, , respectively, with the 
value —1 representing the low level and 1 representing the high level. Letting 
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TABLE 1, 23 + 23 in Regression Form 


Factorial 
Combinations Zo r Zs Te = We Le = Ws Le = Ioly Lr = WTals 


(1) -1 1 1 —1 
—1 —1 7 1 
—1 1 —l 1 
—1 —1 —1 
—l —1 1 
1 —1 —1 
—1 1 —1 
1 1 1 
1 1 —1 
—1 —l -1 
1 -1 —1 
-1 1 


Peet ttt 


the x’s be the elements in matrix X, and defining the observations y as elements 
of a vector Y, we find that 


oOo @ 2 38 4 © © @ 

12 0 0 0 0 0-4 (1) 
0 0 0 O —4 (A) 
12 0 0-4 O (B) 
0 0 O (XY) = (C) 
0 —4 0 O (AB) 
—-4 0 O 12 O (AC) 
0 0 0 0 12 (BC) 

-4 0 0 0 0 0 0 12 (ABC) 


The effects, a vector matrix E, are found by solving the equations given by 
E = (X’X)~'(X’Y). We find, however, that this is readily partitioned into 
four pairs of equations: 


Lael + taal bee ee ae 
ae aluot (ae ee ele 


These are the equations derived in section 4b. The total sum of squares for 
effects, E’ X’ Y also is readily partitioned, as previously indicated. 
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In general when we run a 2" and a 2""* (n = p — q), we have 
a = | 3 le (9) 
(E. 2) +1 3 BR, 
where E, and E, are two factors whose effects, £, and F, , are not clear of one 
another and (£,) and (£,) designate the usual contrasts. The plus sign is used 


when F, and F, are confounded negatively and the minus sign when F, and P, 
are positively confounded. The effects are found by computing 


-1 
Go lel Mie de 
; +1 31 La #,\ l= 3lle@ 

The precisions of the effects may be found as follows: The variance-covariance 
matrix for the contrasts is o”*X’X, so that the variance-covariance matrix for 
pairs of effects is 

2 
V. = (X’X)7'(X’ X)o( XX)! = o°(X’X)' = 


id ort? 


¥1 3 


Thus, the variances of £, and £, are 3 o”/2"**, and their covariance is +o°/2"”. 
The correlation between F, and FP, is therefore p = +3. If the effect were not 
correlated with any other effect, its variance would be o”/3(2"~’). The variance 
of a correlated effect is in the ratio 3 : 4 to the variance o’/2” for only 2” runs. 
Thus, the extra 2"”* runs effect a } reduction in the variance of these effects 
and a 3 reduction in the variance of an uncorrelated effect. 

The 2 d.f. sum of squares is computed by taking 


S.S.for #, and Ff, = B(E,) + BYE, (12) 


If the effect £, is judged real and B£, is not, the sum of squares may be partitioned 
into 1 d.f. sums of squares 


S.S.for 2, = (E,)?/(2" + 2""’) 
S.S. for B, = [F(E,) + 3(E,)]}?/8" + 2””’) 


(13) 


7. A NUMERICAL EXAMPLE 


The data used to illustrate the computations are for a 2°, half of which is 
duplicated. The three factors are designated by the letters A, B, and C. The 
treatment combinations, data, and illustration of the computation of the con- 
trasts are given in Table 2. From the pairs of duplicates we obtain s’ = 2.5425 
with 4 df. 

The analysis of variance, including the decompositions of the two d.f. sums 
of squares for the main effects and two-factor interactions, is given in Table 3. 
It is readily seen that the only effects are due to factors A and B. 

In applying the tests of significance described in Section 4d, we compare 
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TaBLeE 2. Treatment Combinations, Observations and Calculation of Contrasts 


Computation of Contrasts 
Treatment 


Combination Observations Sum II Ill 


—_ 
oo 
bo 


65.4 140.4 291.6 = (J) 
75.0 151.2 26.0 = (A) 
73.2 11.8 14.4 = (B) 
78.0 14.2 --9.6 = (AB) 
—14.8 9.6 10.8 = (C) 
26.6 4.8 2.4 = (AC) 
32.6 41.4 -4.8 = (BC) 
—18.4 -—51.0 -—92.4 = (ABC) 


Bn et 
BrwosynPRre 


do RN Btw 


HO oO 
SRESSRRS 


DNoOoWoONWH 


TaBLE 3. Analysis of Variance 
d.f. Sum of Squares* 


7085. 8800 
(57.7350) 


-o 
Noe Kw e 


- 
— 


~~ 


(22.1400) 


) 


(13.0950) 
~~ 
AB 

ABC 

Error 


2.1600 
10.1700 


1 
1 
2 
1 
1 
1 
4 


Total 


os 
bo 


7191.1800 


* In the decompositions the sums of squares to the left assume that the interaction is zero, 
while the sums of squares to the right assume that the main effect is zero. 


the effects with the standard deviation of an effect, ~/3s°/32 = 0.488. The 
effects are: 


= 2.2875 AC = 0.6750 
= 1.4250 BC = 0.3625 
= 0.7125 ABC = 0.4500 


A 
B 
C 
AB = —0.5625 
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Only A/~/3s"/32 and B/+/3s"/32 exceed 2.78, the value of Student’s distribu- 
tion with 4 d.f. at the 5% level. 


8. EXPERIMENTAL PLANS 


In the prior text we have concerned ourselves with a 2° and a 2°". In this 
section we shall be concerned with designs for more than three variables. 

For four factors we may use a 2* and duplicate half of the treatment combina- 
tions using J + ABCD as generator for the 2*~*. The experimental plan of 24 
runs provides 8 degrees of freedom for error estimation. While in this design 
the three-factor interactions and the four-factor interaction are estimable, we 
might prefer to assume that these interactions are either non-existent or small 
compared to the main effects and two-factor interactions. In this event only 
the two-factor interactions are partially confounded. The general method of 
analysis has been given in Section 6. 

Some difficulties may be encountered with five variables, because of the 
number of runs required by the designs. We may not be willing to make the 
48 runs required for a 2° and a 2°*. A 32 run experiment may be designed by 
completely duplicating a 2°-'. Daniel has suggested a 2°’, half of which is 
duplicated, i.e., a 2°°""*. If we use ] — ABCDE for the 2°* and I — ABCDE + 
AB — CDE for the 2°-'"' then the pairs of effects A, B; C, DE; D, CE; B, CD, 
and 9, AB are partially cnsibiaanibih The: method of analysis ¢ described i in Section 
6 is used. Note that the interactions AC, AD, AE, BC, BD, and BE are con- 
founded only with three-factor interactions. An interesting design given by 
Kempthorne [5] combines a 2°~*, using J — ABCDE as generator, and a 2°”, 
using J — ABC — CDE + ABDE. While this design requires 24 runs, only 4 
degrees of freedom are available for error estimation. 

In the following sub-sections we present experimental designs for the two- 
factor interaction clear (2 fic) fractional factorials with 16 duplicated runs. 
We let p represent the number of 2-level factors and g the degree of fractionation 
of the 2 fic designs. The analysis is given in matrix form for simplicity. 

In the analysis we assume that the factorial combination (1) is always dupli- 
cated, so that any generator or other member of the alias subgroup with an 
odd number of letters is negative. In the analysis we also assume that the inter- 
actions of three or more factors are either non-existent or negligible compared 
with the main effects and two-factor interactions. While this assumption is 
forced with certain of the following designs, the 2’~* is 3 fic, and certain of the 
3-factor interactions are estimable for the other designs. We nevertheless prefer 
to ignore the 3-factor interactions, obtaining a pooled sum of squares by sub- 
traction. 


8a. Six Factors, 32 + 16 = 48 Runs 


The generator for the 2 fic design is ABCDEF; the additional generator 
for the duplicates is — ABC (actually, any three letter generator may be used; the 
use of a two letter generator results in two main effects being partially con- 
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founded). The 48 runs have each of the 6 main effects partially confounded with 
a 2-factor interaction; the nine other 2-factor interactions are clear. 
The general form of the analysis is, for example, 


elles abet “tee 
(BC) —1 3JLBC C 1 34L(BC) 

This form also applies for the pairs of effects B, AC; C, AB; DB, EF; &, DF; P, 
DE. The variance of one of these effects is 3c°/128, a 3 reduction from the 


variance o /32 for the 2 fic design. The other interactions are found by computing, 
for example, 


AD = (AD)/48. (15) 


The variance of these interactions is o”/48, a 3 reduction from the variance in 
the 2 fic design. 
8b. Seven Factors, 64 + 16 = 80 Runs 


The generator for this inefficient design (eight factors require no additional 
runs) is —ABCDEFG; the two additional generators for the 16 duplicates are 
—ABC and —CDE. There are three types of relationship among the effects. 
The groupings and the appropriate analyses are: 


(i) For C, AB, DE, FG: 


1 

5 
7 1 1 (C) 
1 7 —1 —1} (AB) 
1 —l 7 —1] (DE) 
i -] <!] 7IL (FG) 


The variance of each effect is 7 o°/512, a } reduction from the variance o”/64 
for the 2 fic design. 


(i) For A, BC; B, AC; D, CE; E, CD; F, CG; G, CF: 


ola er “ark al = 
(BC) -1 5/.BC BC 1 5JL(BC) 
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(iii) For AD, BE; AE, BD; AF, BG; AG, BF; DF, EG; DG, EF: 
[oo ofe [2] maf] = 8 fA] 
(BE) 1 5ILBE BE —1 5JL (BE) 
For the last two types of relationship the variance of an effect is 5 o°/384, a 
% reduction from the 2 fic design. 


8c. Eight Factors, 64 + 16 = 80 Runs 


The generators for the 2 fic design are -ACDFH and —BDEGH; the two 
additional generators for the 16 duplicates are —ABC and —CDE. There are 
four types of relationships among the effects. The groupings and the appropriate 
analyses are: 

(i) For B, AC; D, CE; F, EG; H, AG: exactly as in equation (17). 
(ii) For AD, BE; AF, EH; BF, DH; BG, CH; CF, DG: exactly as in equation 
(18). 

(iii) For A, BC, GH; C, AB, DE; E, CD, FG; G, AH, EF: 

(A) 5 -1 -1] A A 6 1 1] (A) 

(BC)}=|-1 5 1) BCI, 16-28}BC}=|1 6 -1] BO} «9 

(GH) “4a 1 2a GH} li -1 6{L\@mn 

(iv) For AE, BD, FH; BH, CG, DF: 


(AE) 5 111 AE AE} 6 —1 —11 (4B) 
(BD)|=|1 5 11BD 16-28)BD|=|-1 6 —1]))(BD)} (20) 
(FH) 11 5SiLFH PH -~1 -1  6/.(FH) 


For the first two types of relationship the variance of an effect is 5 o°/384, a 
2 reduction from the 2 fic design. For the last two relationships the variance 
of an effect is 3 o”/324, a + reduction from the 2 fic design. 


_ 8d. Nine Factors, 128 + 16 = 144 Runs 


The generators for this inefficient design (ten and eleven factors require no 
additional runs) are ABDEGH and —ACDEFG_J; the three additional generators 
for the 16 duplicates are -ABC, —EFG, and —ADG. There are only two types 
of relationships among the effects. The groupings and the appropriate analyses 
are: 

(i) For A, BC, DG; B, AC, EH; C, " FJ; D, AG, HJ; E, BH, FG; F, EG, CJ; 
G, AD, EF; H, BE, DJ; J, CF, DH: 


(A) 9 -—1|| A 10 1 = 14 (A) 
(BC) | = 16| —1 111 BC ee BC =!1 10 —1]] (BO| (21) 


(DG) 9 je 1-1 10/(D@) 
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(ii) For AE, CH, DF; AF, BJ, DE; AH, CE, GJ; AJ, BF, GH; BD, CG, 
EJ; BG, CD, FH: 


(AE) 9 1 1) AE AE 10 -1 —1]/(AB) 
(CH)| = 161 9 1\\CH|, 1688\CH}=|-1 10 —1]}(CH)| (22) 
(DF) 1 1 9\LDF DF -1 -1 10/\(~DF/ 


For both relationships the variance of an effect is 5 o°/704, a 1/11 reduction 
from the 2 fic design. 


8e. Ten Factors, 128 + 16 = 144 Runs 


The generators for this relatively inefficient design (eleven factors require 
no additional runs) are -ABEGH, —ACFGK, and ABCDF J; the three addi- 
tional generators for the 16 duplicates are -ADJ, —BCF, and —BDH. There 


are four types of relationship among the effects. The groupings and the appro- 
priate analyses are: 


(i) For A, DH, EF; B, CF, DH; K, CJ, EH: exactly as in equation (21). 
(ii) For AC, BE, DK; AH, BJ, FK: exactly as in equation (22). 
(iii) For C, BF, GH, JK; D, AJ, BH, FG; E, AF, GJ, HK; F, AE, BC, DG; G, 
CH, DF, EJ; H, BD, CG, EK; J, AD, CK, EG: 


(C) § =f =] «J 
(BF) 16 —1 9 1 1 
(GH) —_ -£ § 
(JK) oi 1 
Cc il © 
16-96 BF ied na, 
GH —1 || GH) 
IK 14L(JK) 
(iv) For AB, CE, GK, HJ; AG, BK, DE, FJ; AK, BG, CD, FH: 
(AB) 9 11 114B 
(CE)|_ il 9 1 CE 
(GK) 119 116K 
(HJ) ee a 
ll -1 
—]1 11 
-1l1 -1 
at «il 
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For the first two types of relationship the variance of an effect is 5 o°/704, 
a 1/11 reduction from the 2 fic design. For the last two types of relationship 
the variance of an effect is 11 o°/1536, a 1/12 reduction from the 2 fic design. 


8f. Eleven Factors, 128 + 16 = 144 Runs 


If we wish the variances of the effects to be nearly equal as possible, one of 
the best 2 fic designs has generators —ABEFL, —BCFGK, ADEFGH, and 
BCDEF.J, the three additional generators for the 16 duplicates being —ABC, 
—CDE, and —EFG. (Another design with generators -ABC, —CDE, —EFG, 
—AGH, —AEJ, —CJK, and —GJL for 16 runs is augmentable into a 2 fic 
design with generators —ACDFH, —BDEGH, —BCHJL, and —DFGJK. 
In the 144 run design each of the main effects A, FE, and J are partially confounded 
with five 2-factor interactions, so that there are three “chains” of length six, 
and there are also twelve chains of length four. As we shall see below, the “‘best”’ 
design has six chains of length five and nine chains of length four.) This design 
has three types of relationships. The groupings and the appropriate analyses are: 

(i) For A, BC, DH, EK; F, BJ, DL, EG; G, CL, EF, HJ; J, BF, GH, KL; 
L, CG, DF, JK: exactly as in equation (23). 

(ii) For AF, CJ, GK, HL; AG, BL, DJ, FK; AJ, CF, DG, EL; AL, BG, 
EJ, FH: exactly as in equation (24). 

(iii) For B, AC, DK, EH, FJ; C, AB, DE, GL, HK; D, AH, BK, CE, FL; 

E, AK, BH, CD, FG; H, AD, BE, CK, GJ; K, AE, BD, CH, JL: 
(B) @ «i +] =3§ - 
(AC) -l1 9 1 1 
(DK) | = 16| -1 1 9 i 
(EH) ~—t & €¢ * 
(FJ) 1 1 1 
B aa a ee 
AC a «tf ~) 6 


16-104| DK a. 


EH ce — 
FJ ) af at a ae 


For the first two types of relationship the variance of an effect is 11 o°/1536, 
a 1/12 reduction from the 2 fic design. For the last type of relationship the 
variance of an effect is 3 o”/416, a 1/13 reduction from the 2 fic design. 


8g. Larger Designs 


The 16 duplicate runs effected a } and 3 reduction of the variances of the 
effects with six factors, 3 and + reductions with eight factors, and 1/12 and 1/13 
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reductions with eleven factors. Since the 2 fic designs for 12 to 15 factors require 
256 runs, we would expect even less reduction in the variances. The obvious 
suggestion for more than eleven factors, therefore, is to select a balanced subset 
of the factorial combinations to be duplicated. The correlations among the 
effects should be relatively small, so that averaging the pairs of duplicates 
should not appreciably affect the analysis. 
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Condensed Calculations for 


Evolutionary Operation Programs 


u. E. P. Box ann J. S. HUNTER 
Statistical Techniques Research Group 


Princeton University 


Evolutionary Operation is a method for operating plant processes. Since the method 
is applied by plant personnel themselves as a continuing normal part of process 
operation it is desirable to reduce the calculations to a simple routine matter. This 
paper describes calculation procedures which have been found useful in these circum- 
stances. 


INTRODUCTION 


In 1957 a method for running production processes called Evolutionary 
Operation was described in a paper in Applied Statistics [1]. Since that time the 
method has been adopted extensively in the chemical industry, see for example 
[2], [8]. 


In setting up an Evolutionary Operation program the basic elements are: 


1. provision to introduce a routine of systematic small changes in the levels 
at which process variables are held; 

2. provision to feed back the results derived from making these small changes 
to the operating supervisor responsible for the process. 


A further element which is found to be essential, and which is discussed in 
more detail in [1], is the Evolutionary Operation Committee, an organization 
which continually reviews the results and suggests new action to be introduced 
in later phases. 

The feedback of results is actually accomplished by means of an “Information 
Board”. This may take the form of special blackboard or model placed in the 
supervisor’s office and/or on the factory floor on which an analysis of the data is 
displayed and continuously kept up to date. 

Given the circumstances in which Evolutionary Operation is applied, it is 
essential that the methods for calculating the entries on the Information Board 
be made as simple as possible. Actual experience with a number of Evolutionary 
Operation programs has led the authors to favor a quick and systematic way 
of making the calculations using standard calculation work sheets which entail 
little more than addition and subtraction. Of course, Evolutionary Operation 
depends for its success not only on ensuring that the calculations are easily 
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performed and are likely to be correct, but much more importantly, on a wise 
choice of variables to be studied and on an intelligent appraisal of the results of 
the calculations. There is no absolute way to ensure such things, although 
words of advice and exhortation may help. However, reducing the calculations 
to a simple routine makes more time available for the mental activity needed 
for these more important matters. 

In Evolutionary Operation programs a single performance of a complete 
set of operating conditions is called a cycle and the repeated running through 
of a cycle of operating conditions is called a phase. A new phase of Evolutionary 
Operation begins when new conditions are explored involving different levels 
of the same variables, or different variables. 

In practice several responses (yield, cost, level of impurity, etc.,) would 
usually be portrayed on an information board. In this paper however only a 
single response, the yield (y), is followed in the calculation procedures. For the 
other responses similar calculations would be carried through on separate work- 
sheets. 

In the examples which follow, it is supposed that only a single batch is run at 
each of the sets of operating conditions within a cycle. However, if for some 
reason several batches (runs) are made before changing to a different set of 
operating conditions, the average value of the observations found for each set 
of conditions is treated as a single observation. An Evolutionary Operation 
program may be run on a continuous process by allowing a suitable period of 
time to elapse for the process to settle down after each change in operating 
conditions and then taking one or more observations before instituting a new 
change. The values so obtained are treated exactly like observations from a 
batch process. 

The calculation procedure, which uses the range to provide an estimate of the 
standard deviation and automatically removes block differences, is set out in 
detail here only for two and three variable schemes using two level factorial 
designs with center points. The procedure can however be adapted to any 
number of variables, and to any kind of experimental design. 

In practice, it is often the case that the standard deviation is very similar 
from one phase to the next and, when starting a new phase in which new factors 
are introduced, a good and probably relevant estimate of the standard deviation 
is available from previous phases.* In the examples described here, the previous 
estimate of o is employed in the first cycle or two of each new phase. This 
estimate might not be completely appropriate for the new phase, but it will 
usually be based on a fair number of degrees of freedom so that its sampling 
variation may be ignored and the “limits of error’, +2 (standard error), will 


* At the very beginning of an Evolutionary Operation program, an estimate of the standard 
deviation would often be available from past plant records. However, it is frequently found 
that one consequence of introducing Evolutionary Operation is to produce more careful 
operation of the process and a smaller value of ¢. The use of such a prior estimate is therefore 
often likely to lead to conservative action, 
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provide as reasonable a guide as may be had at this time to the interpretation 
of the results. After two or three cycles a sufficiently reliable estimate of o 
known to be appropriate to the new phase will be obtained, and by this time it 
will again make little difference if the sampling variation in this new estimate 
is ignored. Once again, therefore, the limits of error, +2 (standard error) are 
used. 


CALCULATION PROCEDURE FOR Two VARIABLE PROGRAM 


Suppose that the present phase of an Evolutionary Operation program is 
concerned with studying the effects of time and temperature. The set of k = 5 
operating conditions which together form a single cycle of the Evolutionary 


TEMPERATURE °C 


50 60 70 
TIME (minutes) 


Figure 1 


Operation program are shown in Figure 1. Suppose that the observations listed 
below are obtained in the first four cycles 


Individual Batch Yieids 
Conditions 1 


Cycle I 
Cycle II 
Cycle III 
Cycle IV 


The appearance of the Information Board at the end of the fourth cycle is shown 
in Figure 2. 

The entries shown on the Information Board are the number of the present 
phase, the number of the last cycle completed, the average response observed 
at each set of conditions up to the conclusion of the last cycle, the “error limits” 
for these averages, the “effects” with their associated error limits, the standard 
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INFORMATION Boarp 


_ 
—a—=~osus0s~ouooS3500000wS9S@a{H*VvOo0SSss Oa OOOOeaea+$=$>q00>0000ODDO0" 


Last Cycle Completed 4 


% Yield Other Responses - - - 


Requirement Maximize 


61.9 64.6 
62.2 


63.4 65.7 
_—_—TIME—— 


Error Limits for Averages 


Time 

Effects with Temp 

95% Error ocr 
Limits Change in Mean 


Standard Deviation 


Prior Estimate o 
Figure 2 


deviation calculated from the data of this particular phase, and the prior estimate 
of the standard deviation available from previous work. The “effects” are the 
usual main effects and two factor interactions appropriate for a two level factorial 
design, together with the “change in mean” effect. This latter quantity measures 
the difference between the average yield at the center conditions and the average 
yield at all the conditions being run in this phase. Its use in assessing non- 
linearity and “cost’’ of running the program is described in [1]. 

The entries on the Information Board are changed at the end of each cycle 
as the data become available. The calculations to obtain these entries are most 
conveniently performed on a worksheet such as that shown in Figure 3, using 
a short-cut approximate calculation technique. The worksheets for the present 
example, completed after one, two and three cycles are shown in Figures 3-5. 

The completed work sheet after the first cycle is shown in Figure 3. Since no 
previous observations are available lines (i), (ii), and (iv) on this worksheet 
are left blank. The observations from each set of operating conditions (1), (2), 
(3), (4) and (5) are recorded in line (iii) labelled ‘“New Observations’. In lines 
(v) and (vi), labelled ““New Sums” and ‘“‘New Averages”, are recorded the sum 


Fiacure 3 
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and average of all the observations recorded at each of the operating conditions 
up to and including the cycle considered. In the present case where we are review- 
ing the situation at the end of the first cycle, these sums and averages are simply 
the observations themselves. The standard deviation o cannot be estimated 
from the data of this phase since only one observation is available at each of 
the operating conditions. When, as in this example, a prior estimate of the 
standard deviation is available, this can be employed in calculating provisional 
error limits as indicated in the lower right hand portion of this work sheet. 
At cycle 2 we begin to obtain information from the data themselves concerning 
the value of o relevant to this phase. The steps in filling in the worksheet for 
cycle 2, illustrated in Figure 4, and all subsequent cycles are as follows: 


Step 1: The data from lines (v) and (vi) on the calculation worksheet of the 
previous cycle are copied in lines (i) and (ii) of the worksheet for the 
present cycle. 

Step 2: The observations recorded during the present cycle are listed in line (iii). 

Step 3: The differences between the quantities listed on lines (ii) and (iii) 
are written in line (iv). These are the differences with appropriate 
signs attached between the new observations and the averages over 
previous cycles taken at each set of operating conditions. 

The new sum and new average for each set of conditions are recorded 
on lines (v) and (vi). 

The effects are calculated as given on the worksheet using the new 
averages computed at each of the sets of operating conditions. Thus 
the main effect for the variable ‘‘time” at the end of cycle 2 is given by: 
Time effect = 3(9s + G1 — J2 — Ys) = 2(64.4 + 67.4 — 64.3 — 60.9) = 
3.3. 


After the first cycle it is possible to obtain an estimate of the standard deviation 
of the individual observations (and hence of the error limits for the averages 
and the effects) from the data themselves. These calculations are performed 
in the space provided on the right hand side of the calculation worksheet. 


Step 6: The largest and smallest differences recorded in line (iv) are underlined. 
The range of the differences is the difference between these underlined 
values. This is entered at the right hand end of line (iv). 

Step 7: The range in line (iv) is multiplied by the factor f,,, obtained from 
Table A to obtain an estimate s of o contributed by this cycle. This 
result is entered after ““New s”’ in line (iii). To obtain factor f,,, enter 
Table A with the appropriate cycle number n and appropriate value 
of k (the number of sets of operating conditions in the program). At 
the end of Cycle 2 we have, for k = 5, fs. = 0.30. The estimate of the 
standard deviation contributed by this cycle is then s = Range X fi. = 
4.6 X 0.30 = 1.38. 


Since at cycle 2 no previous estimate of the standard deviation is available 
from the data, the items ‘Previous Sum s” and “Previous Average s” are blank 
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and the entries for ‘‘New Sum s” and “New Average s” are identical with 
that in line (iii). The estimate of the standard deviation after the second cycle 
the prior estimate of ¢ would normally be used in calcu- 
lating the error limits for the “New Averages’, etc. With k = 5 we would 
usually wait until three cycles had been completed before the value of o obtained 
from the data of the phase would be employed. 
All the entries in the calculation work sheet can be filled at the third cycle. 
Using the data given earlier in Table 1, the remaining items will be given for 
the work sheet shown in Figure -5, appropriate at the conclusion of Cycle 3. 


Step 8: The value ‘New s” obtained from step 7 in line (iii) is added to the 
“Previous Sum s” in line (i) and the result recorded on line (v) as 
“New Sum s’”. This quantity is then divided by (n — 1) to give a 
““New Average s’’ in line (vi). 
The error limits for the k averages and the effects are obtained by 
direct substitution of ““New Average s” for o in the equations in the 
lower right portion of the calculation work sheet. 
The new averages, estimated effects and their error limits are recorded 
on the Information Board. The Information Board at the end of Cycle 4 
is shown in Figure 2. 


CALCULATION PROCEDURE FOR THREE VARIABLE PROGRAM 


A single cycle of a three variable Evolutiqnary Operation program follows 
the scheme illustrated in Figure 6. The full cycle is broken into two blocks of 


five runs each as indicated by the open and filled dots in the diagram. The 
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set of operating conditions corresponding to the center of the array is run in 
each block to give a total of ten runs for the full cycle. Operating conditions 
1, 2, 3, 4 and 5 comprise Block I, while 6, 7, 8, 9 and 10 comprise Block II. 

Suppose the observations listed below have been obtained at the end of the 
third cycle. 


BLOCK I BLOCK II 


Conditions 12 3 4 $ a a 


Cycle 1 78 82 63 81 88 85 79 75 78 
2 82 76 79 95 76 69 77 81 70 
3 65 82 68 85 79 86 96 66 83 


The calculation procedure is conducted in two stages: 


(i) the calculation of the main effects, two-factor interactions and “change 
in mean” effect at the conclusion of each cycle, and 
(ii) the calculation of the estimate of the standard deviation at the conclusion 
of each block (twice per cycle). 


It should be noted that there are two calculation work sheets, one for each 
block. Figure 7 shows the situation at the end of the first block, first cycle. 
It is not possible from this half-replicate of the complete design to determine 
the main effects and interactions separately. The quantities FE, , FE, , E; and E£, 
defined on the worksheet may however be calculated. FE, , E, and EF, provide 
estimates of the indicated combinations of main effects and interactions in 
accordance with the well known theory of fractional factorials [4]. EZ, is an 
unbiased estimate of the change in mean effect. On the next calculation work 
sheet are found the corresponding calculations for Block II. From the combined 
information of Blocks I and II, it is now possible to estimate the main effects 
and interactions separately as indicated in the last three lines on the Block II 
worksheet. The error limits for these effects at the conclusion of the first cycle 
are shown on the lower right hand portion of the Block II worksheet using a 
prior estimate of the standard deviation o = 8. 

At Cycle 2 we again begin to obtain information from the data themselves 
concerning the value of o. The steps included on the worksheet parallel very 
closely those used in the two variable scheme and the entries on this and the 
subsequent sheets will be largely self-explanatory. It should be noticed however 
that since there are ten sets of operating conditions (operating conditions 1 
and 6 at the center of the cube are in fact identical but are treated as separate 
since they are performed in different blocks) we must keep tally on ten cumulative 
averages so that, for example, the previous sum for condition 1 will be found on 
the previous worksheet for Block I, and not on the worksheet for Block II. 

From the new averages for the operating conditions 1, 2, 3, 4, 5, we calculate 
new values for E, , E, , E; , and E, which are combined with the new values of 
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E,, E., E, , and E, obtained after completion of Block II, Cycle 2. In this way 
new values for the effects are calculated after each cycle (that is, after each 
pair of blocks). 

In contrast to this, the estimate of the standard deviation is recomputed at 
the conclusion of each block. After the first cycle five new differences are available 
as each new block is formed from which a new estimate of the standard deviation 
can be computed. These estimates of the standard deviation are combined with 
previous values, as shown on the work sheets. 

Using the data recorded earlier for the three variable Evolutionary Operation 
program, the procedures for estimating the effects and their error limits are 
illustrated in the calculation worksheets given in Figures 7 through 12. 


RATIONALE OF COMPUTATION METHOD 


Suppose k sets of conditions are run in an Evolutionary Operation program 
(these might correspond to the sets of conditions in a factorial design or any 
other arrangement whatever). 

At the conclusion of ¢ — 1 cycles a total of k(¢ — 1) observations are available 
which can be classified in a two-way table in which the rows are the ¢ — 1 cycles 
and the columns are the k sets of process conditions. An estimate of the error 
variance having (k — 1)(t — 2) degrees of freedom can be obtained from the 
residual sum of squares S,_; . Consider the data 


Averages for k process conditions 9 , 92, °°: 

at conclusion of (¢ — 1) cycles 
New observations from ¢ cycle a oo ee 
Differences (¥; — y;) Gs his +** ey ss & 


then the following algebraic identity is true: 
AAS) = S, — S,-1 = [(t — 1)/t] * [2d? — (Sd)?/k). 


The A,(S) is the contribution to the residual sum of squares supplied by the 
t® cycle and is distributed quite independently of S,_, . Hence the quantity 
s; = A,(S)/(k — 1) is the unbiased estimate of the variance o” contributed by 
the ¢** cycle alone. In fact an overall estimate of the variance after n cycles is 
of the additive form 


’ = : AAS)/(k — Din — 1) = > s/(n — 1} 


Now as we have already mentioned there is usually no lack of data from which 
the standard deviation can be calculated. There seems little point therefore in 
not using the slightly less efficient estimate of the standard deviation obtained 
from the range of the k differences d,; , d, , --- , d, . In both the two- and three- 
variable schemes discussed above the value of k is five, and for such values, the 
estimate of « obtained from the range is known to be very little less efficient [5], 
and somewhat more robust [6] than the estimate based on the sums of squares. 
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To see how the range is used we notice that if we write 
2; = [(t — 1)/t}'d; 


where d; is the 7‘ difference in the ¢** cycle. Then 


Si a » [((é om 1)/t)(d; a a)’ = a (z; si z,)” 


t=1 


An appropriate alternative estimate of the overall standard deviation from 
n cycles would therefore be provided by 


sm Fae” 


where §, was obtained from the range of the z; . However, the range of the z; 
is [(t — 1)/t]'’w, where w, is the range of the differences d; in the ¢* cycle. 
Thus our estimate §, is given by 


& = [(t — 1)/t]'ew, 


where c, is the proper constant (obtained for example from the Biometrika 
Tables for Statisticians [5}) by which the range w, must be multiplied to provide 
(on the assumption of Normality) an unbiased estimate of o from k observations. 

If we write f,.. = ¢,{(n — 1)/n]'” then s, = f,,nW, . The value f,,, for various 
values of k and n are given in Table A. At the conclusion of the n** cycle there 


TaBLeE A 
Values of Constant fin 


Number of k = Number of sets of conditions in block 
cycles = n 3 4 5 


34 
40. 
42 
43 
44 
45 
45 
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are n — 1 estimates. The combined estimate of the standard deviation used 
for the Information Board is simply the average of these. 


We are grateful to J. W. Tukey for valuable suggestions concerning the 


organization of this material, and to Miss Gillian Richardson who checked the 
calculations. 
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NOTICES 


NaTIONAL MEETINGS OF THE AMERICAN SOCIETY FOR QUALITY CONTROL 


The Thirteenth Annual Convention of the American Society for Quality 
Control and the National Quality Control and Production Exposition will be 
held this year on May 25, 26 and 27 in Cleveland, Ohio. The convention head- 
quarters will be the Hotel Cleveland, with the convention program being held in 
the Cleveland Public Hall. The chairman of the meetings is Mr. Wade R. 
Weaver. Seventy-five papers are presently scheduled for the technical sessions. 
The technical sessions will also include a Department of Defense Quality 
Control Symposium under the direction of John J. Riordan, Staff Director, 
Inspection and Quality Control Division, Office of the Secretary of Defense. 
This symposium will include sessions titled ‘Military Progress in Quality 
Control in 1958” with Mr. Riordan as chairman; ‘New Developments in 
Military Food Inspection” under the chairmanship of Col. Richard L. Lewis, 
Chief, Inspection Division, Military Subsistence Agency; and a third session on 
“Quality Control and Reliability of Military Equipment” under the chairman- 
ship of Dr. Max Astrachan, Air Force Institute of Technology, Wright-Patterson 
Air Force Base and Mr. John J. Crowley of the General Dynamics Corporation. 

Complete details of the program of the convention will appear in Industrial 
Quality Control. For further information concerning the meetings, or the Exposi- 
tion, write Mr. Wade R. Weaver, Republic Steel Corporation, Cleveland, 1, Ohio. 


SOUTHERN REGIONAL GRADUATE SUMMER SESSION IN STATISTICS 
AT Nortu CAROLINA STaTE COLLEGE 


The 1959 session of the Southern Regional Graduate Summer Session in 
Statistics will be held at North Carolina State College, Raleigh, from June 8 
to July 17, 1959. North Carolina State College, Virginia Polytechnic Institute, 
University of Florida, and Oklahoma State University have operated this 
continuing program of graduate summer sessions in statistics since 1954. 

The 1959 session, like previous sessions under this program, is intended to 
serve: 1) teachers of introductory statistical courses who want formal training 
in modern statistics; 2) research and professional workers who want intensive 
instruction in basic statistical concepts and modern statistical methodology; 
3) professional statisticians who wish to keep informed about advanced specialized 
theory and methods; 4) prospective candidates for graduate degrees in statistics; 
and 5) graduate students in other fields who desire supporting work in statistics. 

The 1959 session will last six weeks and courses will carry three semester hours 
of credit. Not more than two courses may be taken for credit at any one session. 
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The summer work in statistics may be applied as residence credit at any of the 
cooperating institutions, as well as certain other universities, in partial fulfillment 
of the requirements for a graduate degree. The program may be entered at any 
session, and consecutive courses will follow in successive summers so that it 
would be possible for a student to complete the course work in statistics for a 
Master’s degree in three summers. Students must satisfy the remaining require- 
ments for course work and thesis at the institution where they are to be admitted 
to candidacy. The advanced courses may be accepted as part of the Ph.D. 
program of the participating institutions. 
The courses to be offered in statistics in 1959 at Raleigh are as follows: 


Methods I. Distribution of sample statistics, #, s°, x’, ¢ and F, and their 
uses in estimation and tests of hypotheses; regression; correlation; elementary 
experimental and sample design; analysis of variance including multiple com- 
parison techniques. 

Methods II. Principles of scientific experimentation, covariance, multiple 
regression, factorial experiments, confounding, incomplete blocks, balanced and 
partially balanced. 

Theory I. Probability, probability distribution for discrete and continuous 
variates, expected values and moments, sampling distributions, estimation, 
tests of significance. 

Theory II. Continuation of Theory I, includes sampling distributions of s’, 
t, F; tests of hypotheses; interval estimation; goodness of fit. 

Theory III. A study of the theory of multiple linear regression, estimation 
and tests of hypotheses of parameters in multiple regression, Markoff theorem, 
distribution of quadratic forms, power of tests. 

Stochastic Processes. Random walks; Markov chains in discrete and continuous 
time; Markov processes; diffusion processes; first passage times; renewal theory; 
regenerative theory; stationary processes; spectral and prediction theory; evolu- 
tionary stochastic processes; applications to population growth, stores, queues, 
inventories, etc. 

Methods of Operations Research. Linear Programming, theory of games, 
techniques for analyzing waiting lines and queues; applied probability; recent 
developments, applications of results to specific problems; case studies. 

Theory of Sampling. Basic theory of sampling from a finite population. Con- 
fidence limits and estimation of optimum sample size, comparison of different 
sample designs, methods and probabilities for selection and methods of estima- 
tion, choice of a sampling unit, double sampling, matched samples. 

Advanced Statistical Analysis. General computational methods for linear 
regression, non-orthogonal data, carryover effects, orthogonal polynomials, 
response surfaces, non-linear systems, variance components for orthogonal and 
non-orthogonal data. 


A number of courses will, in addition, be available in the Mathematics Depart- 
ment. Courses to be offered during the June 8-July 17 period of the Summer 
Session in Statistics include: 
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Differential Equations; 

Introduction to Determinants and Matrices; 
Introduction to Applied Mathematics; 
Numerical Analysis; 

Advanced Calculus I; and 

Boundary Value Problems. 


The formal course offerings will be supplemented by seminars and special lectures. 
Requests for application blanks, and for further information should be 
addressed to 
Professor F. E. McVay 
Department of Experimental Statistics 
North Carolina State College 
Raleigh, North Carolina 


NATIONAL SCIENCE FouNDATION SUMMER INSTITUTE OF STATISTICS FOR 
CoLLEGE TEACHERS OF MATHEMATICS 


An announcement has been made by the National Science Foundation of a 
Summer Institute for College Teachers of Mathematics specialized in the subject 
matter of Statistics to be held at Oklahoma State University, June 8 through 
July 31, 1959. The objective of this program is to provide college teachers of 
mathematics, whose primary competence is in some field of mathematics other 
than statistics, an opportunity to obtain basic training in the theory of statistics 
and probability, and to experience the use of statistical methods in scientific 
research. The Institute is designed to provide the motivation and background 
necessary for a teacher of mathematics to perform in a creditable manner in 
teaching statistics, either theory or methods, on an undergraduate level. 

There are three basic courses being offered which will be equivalent to eight 
semester hours of graduate credit. Instruction will be organized into lecture, 
laboratory and seminar periods all conducted in air-conditioned classrooms. 


Course 1: Probability and Distribution Theory 
Course 2: Theory of Statistical Inference 
Course 3: Selected Topics in Statistical Methods. 


Special seminar speakers who have tentatively accepted include 


Irving Burr Purdue University 

Dudley Cowden University of North Carolina 
William Krumbein Northwestern University 

William G. Cochran Harvard University 

Morris Hansen Bureau of Standards 

A. W. Wortham Texas Instruments 

Raymond J. Jessen General Analysis Corporation 

S. M. Free Smith Kline and French Laboratories 
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The Institute is limited to thirty participants. Requests for further information 
should be addressed to: 


Dr. Carl E. Marshall 

Director, Statistical Laboratory 
Oklahoma State University 
Stillwater, Oklahoma. 


GoRDON CONFERENCE ON STATISTICS IN CHEMISTRY AND CHEMICAL 
ENGINEERING 


The Gordon Conference on Statistics in Chemistry and Chemical Engineering 
will be held again this year at the New Hampton School, New Hampton, New 
Hampshire, from the 24th through the 28th of August, 1959. The purpose of 
each Gordon Conference is to “extend the frontiers of science by fostering a 
free and informal exchange of ideas among persons actively interested in the 
subjects under discussion. The purpose of the program is not to review the known 
fields of chemistry but, primarily, to bring experts up to date on the latest 
developments, to analyze the significance of these developments, and to provoke 
suggestions concerning the underlying theories and profitable methods of 
approach for making new progress.” 

The chairman of this year’s conference on Statistics in Chemistry and Chemical 
Engineering is Dr. R. J. DeGray, The Standard Oil Company (Ohio), Cleveland, 
Ohio. Topics and speakers presently scheduled for the conference include: 


Evolutionary Operation F. J. Riordan 
Experimental Designs to Adjust for Time Trends Hubert Hill 
Partially Duplicated Factorials : Otto Dykstra 
Multiple Regression Analysis H. C. Hamaker 
Non-Linear Estimation G. E. P. Box 
Curve Fitting Gabriel Kron 
Random Balance James Dolby 
Fixed, Mixed and Random Model C. A. Bennett 
Change-Over Trials H. L. Lucas 


Those wishing to attend the conference should write to: 


Dr. W. George Parks 
Department of Chemistry 
University of Rhode Island 
Kingston, Rhode Island 


Attendance at the conference is by application only, and the number of attendees 
that can be accepted for the conference is limited. For complete information on 
the 1959 Conferences, see Science, 6 March 1959. 


The paper ‘Measurements Made by Matching with Known Standards’ by W. J. Youden, 
W. S. Connor and N. C. Severo was accepted for publication in this issue but will appear in 
the next issue instead. 
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PREPARATION OF MANUSCRIPTS 


Manuscripts should be submitted to the office of the editor: J. 8. Hunter, 167 
Nassau St., F’rinceton, New Jersey. Each manuscript should be typewritten, 
double spaced, with wide margins at sides, top, and bottom. The original should 

be submitted with two additional copies, on paper that will take corrections. 
Tousah it WAMUGcMGA copes ame omagunie cle eomebady teabae 
Footnotes should be avoided and replaced by remarks in the text, or placed in 
an appendix. Preferably, references in the manuscript should appear as (Jones, 
4. B., 1958), and again later in alphabetical order in a list of references. Al- 
ternatively references may be numbered, e.g. [1], as they appear in the manu- 
script and be listed in this sequence in the list of references. In the reference 
list, each reference should contain, in the order indicated, the name and initials 
of the author followed by those of the co-authors, date of publication, title of 
reference, source, volume number and page. References to books should include 
publisher’s name and location. 

Figures, charts, and diagrams should be professionally drawn on plain white 
paper or tracing cloth in black: India ink twice the size they are to be printed. 
A iull page diagram, in print, measures 7.25 X 4.75 inches. 

As far as possible, formulas should be typewritten and symbols not available 
on a typewriter carefully inserted in ink. Authors are asked to keep in mind the 
typographical difficulties of complicated mathematical formulae. The difference 
between capital and lower-case letters should be clearly shown; care should be 
taken to avoid confusion between such pairs as zero and the letter O, the numeral 
1 and the letter 7, numeral 1 used as superscript and prime (’), alpha and a, kappa 
and k, mu and u, nu and », eta and n, etc. Subscripts or superscripts should be 
clearly below or above the line. Bars above groups of letters (e.g., log z) and 
underlined letters (e.g., z) are difficult to print and should be avoided. Symbols 
are automatically italicized by the printer and should not be underlined on 
manuscripts. Boldface letters may be indicated by underlining with a wavy line 
on the manuscript; boldface subscripts and superscripts are not available. 
Complicated exponentials should be represented. with the symbol exp particu- 
larly when appearing in the text, that is, 


exp [(@’ +- &°)] should be used in place of e‘*"**”*”*, 


In writing square roots the fractional exponent is preferable to the radical sign. 
Fractions in the body of the text (and when possible in displayed expressions) 
and fractions occurring in the numerators or denominators of fractions are 
preferably written with the solidus; thus 


a+b 
(a + b)/(c + d) rather than ea 


Authors will ordinarily receive only galley proofs. Fifty offprints without 
covers will be furnished free. Costs for additional reprints and covers can be 
furnished on request. 
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